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A  Navier-Stokes  problem  solver,  developed  by  L.  N.  Sankar.  is  modified  to  provide 
dynamic,  interactive  graphical  presentations  of  predicted  flow  field  solutions  about  a 
NACA-0012  airfoil  section,  oscillating  in  pitch,  in  order  to  demonstrate  the  capabilities 
of  dynamic  graphics  applications  in  the  study  of  complex,  unsteady  flows.  Flow  field 
solutions  in  the  form  of  pressure  coefficient  and  stream  function  contour  plots  about  an 
airfoil  experiencing  dynamic  stall  are  plotted  utilizing  an  IRIS  3000-series  workstation 
and  Graphical  Animation  System  (GAS)  software,  developed  by  Sterling  Software  for 
NASA.  These  full  cycle  solutions,  in  conjunction  with  dynamic  surface  pressure  dis¬ 
tribution  plots  and  integrated  lift,  pitching  moment  and  drag  coefficient  data,  are  com¬ 
pared  to  existing  experimental  data  in  order  to  provide  an  indication  of  the  validity  of 
the  code's  far-field  solution.  Full  procedural  documentation  is  maintained  in  order  to 
provide  an  efficient  analysis  tool  for  use  in  future  oscillating  airfoil  studies  planned  by 
the  NASA-Ames  Fluid  Mechanics  Laboratory  and  the  Naval  Postgraduate  School  De¬ 
partment  of  Aeronautics  and  Astronautics. 
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I.  INTRODUCTION 

It  has  been  recognized  for  some  time  [Ref.  1]  that  an  airfoil  oscillating  m  pitch  to 
angles  of  attack  greater  than  the  static  stall  angle  will  surpass  the  traditional  stud  umicr 


and  generate  normal  forces  which 


:ed  those  attainable  in  i he  '.aiic  case.  I  h.s  dy¬ 


namic  stall  phenomenon  is  attributed  to  the  aft  movement  of  a  strong  shed  voicex  .  long 
the  upper  surface  of  the  ah  Toil,  carrying  with  it  an  induced  velocity  field  whi.h  radically 
and  dynamically  changes  chordwise  pressure  distiibutions.  In  geneial.  an  accurate  in¬ 
terpretation  of  the  dynamic  stall  mechanism  will  significantly  impact  a  variety  cl  appli¬ 
cations,  all  of  which  involve  dynamic  lifting  surface  motions  and  unsteady  few 
separation.  Originally,  dynamic  stall  analysis  efforts  were  diiected  toward  helicopter 
aerodynamics,  where  sharp  increases  in  oscillatory  torsional  loading  and  thus,  blade 
stress,  can  reduce  the  fatigue  life  of  rotor  mechanical  components  and  vortex-induced 
aerodynamic  loading  can  generate  adversely  phased  pitching  moments,  resulting  m  stall 
flutter  [Ref.  2|.  Much  current  interest  concerns  the  feasibility  of  exploiting  dynamic  stall 
forces  for  effective,  sustained  maneuvering  in  the  high  angle  of  attack,  "post-sta'l"  flight 
regime,  or  supermancuverability.  for  next-geneiation  fighter  and  attack  aircraft  [Ref. 
Historically,  attempts  to  analyze  such  complex,  unsteady  behavior  relied  heavily  on 


empirical  data,  obtained  from  often  laborious  and  time-consuming  test?.  With  the  ad¬ 
vent  of  the  supercomputer,  however,  computation  of  actual  viscous  flew  held?  about 
moderately  complex  computational  models  can  now  be  numerically  achieved  in  a  matter 
of  minutes  through  solution  of  the  Reynolds-averaged  Navier-Stokes  equations  [P.ef. 
4).  These  solutions,  in  and  of  themselves,  however,  are  not  sufficient  to  promote  insight 
into  the  mechanics  and  physics  involved  in  such  flows.  This  additional  requirement,  lor 
effective  visual  portrayal  of  the  How  field,  is  satisfied  by  application  of  high  performance 
interactive  computer  graphics  workstations  and  associated  software. 

The  long  range  goal  of  the  Compute  tional  Fluid  Dynamics  held  is  the  development 
of  a  thoroughly  verified  computer  code  for  unsteady  aerodynamics  which,  among  a 
myriad  of  other  applications,  will  provide  future  aircraft  designers  with  the  opportunity 
to  derive  full  advantage  from  application  of  the  dynamic  stall  phenomenon.  To  this  end. 
the  Tluid  Mechanics  Laboratory  (FML)  of  the  NASA -Ames  Research  Center  (ARC), 
in  conjunction  with  the  Naval  Postgraduate  School  Department  of  Aeronautics  end 


1 


Astronautics  has  planned  an  assortment  of  parallel  and  complementary  oscillating  airfoil 
windtunnel  experiments  and  computer  simulations. 

The  current  study  utilizes  existing  dynamic  graphics  packages  to  generate  real-time 
projections  of  How  fields  about  a  NACA-O012  airfoil  oscillating  in  pitch  and  experienc¬ 
ing  deep  dynamic  stall,  as  provided  by  a  Navier-Stokes  solver  developed  by  L.  N.  Sankar 
[Refs.  5,h).  A  modified  version  of  the  Sankar  code  was  submitted  via  the  FMI.  front  end 
VAX,  to  the  NASA-ARC  Cray  X-MP  48  computer,  which  output  flow  field  solutions 
at  specified  intervals  throughout  the  oscillatory  cycle.  From  this  data,  graphics  files  were 
generated  which,  in  turn,  were  submitted  to  the  Graphics  Animation  System  (GAS) 
software  as  developed  by  Sterling  Software  under  contract  to  NASA-ARC,  and  run  on 
an  IRIS  3000-series  graphics  workstation.  (Final  output,  for  demonstrative  purposes, 
was  then  transferred  to  video  tape.)  Thus,  the  results  of  the  studs  :re  two-fold.  First, 
procedures  by  which  interactive  computer  graphics  could  be  efficic.  (in  terms  of  both 
computer-time  and  man-hours)  and  effectively  (in  terms  of  information  display)  incor¬ 
porated  as  an  analysis  and  verification  tool  for  future  FML  studies,  were  developed, 
tested  and  refined.  Secondly,  the  resultant  graphics  were  utilized  as  a  tool  for  the  on¬ 
going  verification  of  the  Sankar  code. 


\rm  '^. 


II.  DESCRIPTION  OF  THE  SANKAR  CODE 


Simulation  of  complex  phenomena  occurring  in  real  fluid  flows  requires  accurate 
solutions  to  the  full  Navier-Stok.es  equations.  The  Sankar  Naviei -Stokes  solver.  dev el¬ 
oped  for  Bfade-Vortex  Interaction  (BYh  studies,  solves  the  two-dimensional,  unsteady, 
compressible  fiulcr  and  Navier-Stokes  equations  in  strong  conservation  form,  utilizing 
an  alternating  direction,  implicit  method  as  the  time  marching  algorithm.  A  body-fitted 
C-grid.  with  clustering  in  the  normal  diiection  is  utilized  to  discretize  the  How  field. 
Turbulent  shear  stresses  are  simulated  with  a  two-layer  algebraic  eddy-viscosity  model, 
with  modifications  as  described  later.  The  code  may  thus  be  utilized  for  solution  of 
steady  or  unsteady,  inviscid  or  viscous  and  laminar  or  turbulent  flows. 


A.  GOVERNING  EQUATIONS 

In  Cartesian  coordinates,  the  2-D,  unsteady,  compressible  Navier-Stok.es  equations 
in  stronc  conservative  form  mav  be  written  as 


d,q  -f  SXE  4-  S  F  =  Re  (SXR  +  S  S) 


(2.11 


where  q,  £,  F,  R  and  S  are  four  element  vectors  with  entries  corresponding,  in  order,  to 
the  equations  of  continuity,  momentum  (x-  and  y-)  and  energy.  If  non-dimensionaiized 
such  that  all  values  of  length,  velocity  (u  and  v),  density  (  p  )  and  total  energy  per  unit 
volume  (e)  are  normalized  with  respect  to  section  chord  length  (c),  free  stream  speed  of 
sound  (  cix ),  free  stream  density  (  px)  and  p„ux,  respectively,  these  vectors  may  be  writ¬ 
ten  as 
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The  Reynolds  number  (Re),  pressure  (p),  speed  of  sound  (a)  and  stress  terms  (  r,.)  are 
defined  as  follows. 
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(2.3 

ly  —  1  l[e  —  0.5/.)(u*  4-  v‘)J 

(2.4) 

;■(’,'  -  1  )[c/p  -  0.5(u'  +  V2)] 

(2.5, 

Txx  =  ('•  +  2p)ux  +  /  l\. 

(2.6; 

r<y  =  ."("v  +  vA) 

Tyy  =  0.  +  2fl)vy  +  >UX 

=  tiTxx  +  vrxy  +  kPr~\y  -  \  )SXU2 
54  =  urxy  4-  v-yy  +  kPr~\y  - 1  )6ya2 

Assuming  Stokes'  hypothesis  to  be  valid,  the  constant  /  is  defined  as  -2  3  m  .  Since  only 
airflow  is  considered,  the  ratio  of  specific  heats  (  y  )  is  defined  as  1.4  and  the  Prandtl 
number  (Pr),  as  one. 

Neglecting  viscous  terms  results  in  the  Euler  equations,  as  the  right-hand  side  of 
equation  (2.1)  goes  to  zero.  Use  of  the  strong  conservative  form  (i.e.,  the  continuity 
entry  in  E,  6(pu)idx,  is  solved  in  its  present  form,  vice  the  more  mathematically  correct 
form  of  pdulo.x  +  udpldx)  allows  the  identical  conservation  of  physical  quantities  (mass, 
momentum  and  energy)  when  finite  difference  schemes  are  applied. 

B.  TRANSFORMED  GOVERNING  EQUATIONS 

1  he  flow  field  region  must  be  discretized  into  a  transformed,  finite  difference  mesh, 
or  computational  plane,  in  order  to  allow  numerical  solution  of  the  governing  equations. 
The  most  efficient  grids  are  rectangular  in  shape  with  regular  spacing  or  spacing  gradi¬ 
ents.  They  must,  however,  correspond  to  a  flow  field  grid  which  provides  high  resolution 
(minimal  spacing)  in  regions  where  gradients  are  large,  particularly  in  the  boundary  layer 
and  about  the  leading  edge.  In  response  to  the  coordinate  transformation  involved  in 
the  grid  generation  process,  the  governing  equations,  too,  must  be  transformed.  By  de¬ 
fining  s  and  as  functions  of  the  cartesian  coordinates  (time  is  not  transformed),  this  is 
simply  a  2-D  mapping  procedure  accomplished  by  application  of  the  transformation 
Jacobian, 
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By  setting 


h(.t\t)  oxl dxjdq 

^U>  V)  Syj S’  SyjSij 


<i( ») 
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the  governing  equation  becomes 

(5. <7  4-  d-£  +  S^F  =  Re  '(<)-/?  +  b);S) 

where 

q  =  q!J 

E-(ti  +  ZxE  +  ZyF)IJ 
F={ri,q  +  >/,£  +  >iyF)IJ 
R  =  +  C: yS)IJ 

S  ~  (qxR  +  *tyS)U 

and 

tjuc  =  (>•  +  2m)(Cxw;  +  */*“,)  +  'i-(Cy‘{  +  Vl 
Txy  =  m[(7mC  +  VyUr)  +  (Cxv;  + 


Tw  =  ('-  +  2^yv;  +  Vl  +  ^(CxMC  + 

=  wrxx  +  VTxy  +  kPr~\y  -  l)(Cxdca2  +  qxSna2) 

SA  =  urxy  +  VTyy  +  kPr~\y  -  1  KC^-a2  +  v/^a2) 


C.  GRID  GENERATION 

Grid  generation  for  discretization  of  the  How  Ik lu  region  may  be  aocompli'hed  by 
conformal  mapping,  algebraic  methods  nr  numeiitallv  solving  a  set  of  partial  Jdieromi.i: 
equations.  The  original  Sankar  code  utilizes  either  of  the  latter  two  methods,  uhile  the 
version  currently  vectorized  for  use  in  this  study  utilizes  only  the  the  algehiaie  method, 
which  results  in  a  body-litted.  sheared  parabolic  coordinate  system  or  (.'-grid.  I  nh.  a'ic  n 
ol  such  body-litted  grids  allows  synchronous  rotation  of  the  entire  g:ul  and  .uifoil  sec¬ 
tion  for  dynamic  cases,  using  simple  trigonometric  relations.  This  coordinate  system 
satisfies  the  general  gi id  requirements  for  smoothness  throughout  and  line  -paune  in 
regions  where  high  gradients  exist,  such  as  the  boundary  layer  or  leading  edge. 

Normalized  geometric  airfoil  shape  data,  in  Cartesian  coordinates,  is  input  m  umle 
format  to  the  code,  which  utilizes  an  interpolative  procedure  to  compute  additional 
points  and  smooth  the  surface.  I  he  trailing  edge  region  is  modelled  as  a  vortex  sheet 
shape,  or  "cut",  which  smoothly  leaves  the  airfoil,  tangent  to  the  mean  camber  line  at 
the  trailing  edge.  Algebraic  manipulation  of  the  section  surface  and  cut  allo  ws  grid 
generation  in  the  transformed  I  -  1/  plane.  Figure  1  shows  the  resultant  grid  when 
mapped  back  to  Cartesian  coordinates. 

Uniform  spacing  in  the  transformed  plane's  f  -direction  results  111  line  spacing  at  the 
leading  edge  in  the  real  plane,  but  coarse  spacing  in  the  nailing  edge  legion,  due  to 
uniform  increments  across  the  axis  (as  opposed  to  the  standard  C-grid  which  results  in 
a  region  of  finer  spacing  at  the  trailing  edge).  >1  -spacing  in  both  planes  increases,  ap¬ 
proximately  exponentially,  in  directions  normal  to  the  aii foil  surface,  with  the  magnitude 
of  initial  spacing  being  user  defined.  Though  easy  to  generate,  this  grid  s  effectiveness 
in  the  analysis  of  certain  flow  cases  has  proven  somewhat  limited,  due  to  its  coarseness 
in  the  trailing  edge  region  [Ref.  6,  p.44[. 

D.  APPLICATION  OF  BOUNDARY  CONDITIONS 

Consideration  of  an  airfoil  impulsively  started  from  rest  in  a  lluid  with  uniiorm 
properties  throughout,  provides  the  initial  conditions  for  solution  of  the  parabolic  Liulcr 
and  Navier-Stokes  equations  (Eq.  2.9)  in  the  computational  piano.  In  addition,  bound¬ 
ary  conditions  and  artificial  dissipation  terms  are  required  in  order  to  achieve  accuiate 
solutions. 

Three  boundaries  exist  which  are  of  interest,  the  surface,  the  far-held  boundary  lgiid 
limit)  and  the  trailing  edge  cut.  Boundary  conditions  at  the  surface  are  driven  by  the 
110-slip  condition  which,  in  response  to  viscous  elfects,  requires  the  (laid  velocity  at  the 


boundary  to  equal  the  boundary  velocity.  Thus,  in  this  reference  frame,  u  and  v  .tie  sec 
to  zero  on  the  solid  surface.  Since  the  surface  is  considered  adiabatic,  both  <v  a-/  and 
Jt’/i)'  are  set  to  zero  and.  likewise,  pressure  distributions  aie  determined  from  op  os;  and 
Cpj o'  equal  to  zero  at  the  solid  boundary.  (This  is  equivalent  to  neglectme  scree;  con¬ 
tributions  in  the  momentum  equations,  which  is  acceptable  when  dealing  with  high 
Reynolds  numbers  flows.)  Boundary  conditions  at  the  cut  are  driven  by  the  necessity, 
to  avoid  discontinuities  in  the  "continuous"  poition  of  the  How  field.  Since  the  end  is 
extremely  dense  in  the  1/  -direction  in  this  region,  averages  of  the  values  of'  the  two 
nearest  interior  points  are  assigned  to  points  along  the  cut,  allowing  a  smooth  transition. 
Since  the  grid  cannot  economically  be  generated  large  enough  to  avoid  disturbance  ve¬ 
locities  at  the  fur-lield  boundary,  boundary  conditions  must  account  for  the  presence  of 
the  airfoil.  Linear  small  disturbance  theory  is  applied  to  determine  perturbation  veloci¬ 
ties  at  points  along  the  boundary,  which  are  then  added  to  free  stream  conditions. 
Downstream  boundaries  are  treated  such  that  entropy  changes  can  be  confected  out  of 
the  computational  domain  [Ref.  6.  p.  29j,  allowing  shocks  and  boundary-layer  generated 
vorticity  to  pass  through  the  grid. 

It  has  been  found  that  spatial  derivatives  are  sensitive  to  the  decoupling  of  odd  and 
even  points  which  necessarily  occurs  in  central  difference  schemes.  This  results  in  the 
generation  of  high  frequency  errors  in  regions  of  large  pressure  gradients,  sm.li  as  are 
present  in  shocks  or  about  stagnation  points.  Due  to  the  high  Reynolds  numbers  en¬ 
countered  in  these  flows,  dissipation  provided  by  the  viscous  terms  arc  not  enough  to 
eliminate  the  errors,  necessitating  the  addition  of  two  artificial  dissipation  terms,  em¬ 
bedded  in  the  numerical  schemes.  An  explicit  artificial  viscosity  term  is  input  by  the 
user,  with  a  proportional  implicit  term  assigned  by  the  code. 

E.  TURBULENCE  MODELLING 

Since  during  the  derivation  of  the  Nav  ier-Stokes  equations,  no  assumptions  are 
made  regarding  flow  type,  these  equations  are  instantaneously  valid  for  both  laminar 
and  turbulent  flows.  The  large  range  of  time  and  spatial  scales  encountered  in  turbulent 
flows,  however,  makes  solution  of  the  instantaneous  Navier-Stokes  equations  viituully 
impossible,  at  this  time.  As  a  result,  the  equations  are  time-  or  ensemble-averaged.  Use 
of  such  equations,  in  response  to  turbulence,  results  in  additional  turbulent  shear  stress 
terms  or  the  Reynolds  stresses  (named  for  Oswalde  Reynolds,  who  initiated  turbulence 
studies  in  the  ISSO's),  which  effectively  are  increases  in  shear  stresses  due  to  turbulent 
motion.  The  difficulties  associated  with  the  existence  of  additional  unknowns  without 
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any  additional  equations  is  described  as  the  closure  problem  and  is  dealt  with  through 
the  use  of  turbulence  models.  These  are  based  on  empirical  knowledge  of  turbulent 
flows  and  thus,  provided  solutions  will  be  approximate  and  require  some  confirmation 
from  experiment  [Ref.  7|. 

The  most  common  turbulence  modelling  method  involves  mixing  lengths,  as  devel¬ 
oped  by  Prandtl.  Assuming  that  turbulent  fluctuations  are  essentially  the  result  of  ve¬ 
locity  perturbations  between  adjacent  streamlines,  and  that  all  fluctuating  velocity 
components  at  a  given  point  are  of  the  same  magnitude,  it  can  be  shown  that  turbulent 
or  eddy  viscosity  (  /ur  )  is  proportional  to  the  magnitude  of  the  local  vorticity  (  u  )  [Ref. 
8,  p.  3S7],  In  Prandrls  mixing  length  model,  the  proportional  term  is  the  square  of  a 
characteristic  length,  related  to  fluid  turbulence  intensity.  Prandtl  suggested  this  char¬ 
acteristic  length  (  /  )  be  treated  as 


l  =  ky 


where  y  is  the  normal  distance  from  the  fluid  boundary  and  k  is  empirically  obtained. 
Thus,  the  key  to  obtaining  accurate  solutions  when  utilizing  models  of  this  type  lies  in 
the  mixing  length  expression. 

The  Baldwin-Lomax  two-layer  eddy  viscosity  model,  based  on  Cebeci's  two-layer 
model,  is  presently  used  in  the  Sankar  code.  This  model  divides  the  boundary  layer  into 
two  regions,  an  inner  layer  and  an  outer  layer,  with  separate  methods  for  determining 
eddy  viscosity  used  in  each.  The  boundary  for  the  two  layers  is  defined  as  that  point 
where  eddy  viscosities  produced  by  the  two  methods  match.  1  he  inner  layer  utilizes  a 
mixing  length  model  where  eddy  viscosity  starts  from  zero  at  the  wall  and  is  defined  by 
the  following. 

(Hr)inner  =  pl2  M  (2.12) 

The  mixing  length  term  is  defined  by 

/  =  k}D  (2.13) 

where y  is  the  normal  distance  from  the  wall,  k  is  the  Von  Karman  constant  and  D  is  the 

Van  Driest  damping  function,  given  by 

D  =  [l  -exp(-/M+)]  (2.14) 

where  A~  is  an  empirical  constant.  The  outer  layer  eddy  viscosity  model  is  defined  by 
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where  K  and  Ccp  are  constants.  The  K  motf  intermittency  function,  given  by 

Wf)  =  [1  +  5-HCK,„by!ym3Xfr'  l-l 6) 

ensures  that  eddy  viscosity  approaches  zero  as  the  edge  of  the  boundary  layer  is  ap¬ 
proached  and  the  flow  assumes  external  characteristics.  CKM  is  a  constant.  l\ak,  is  a 
function  drilled  the  following  relation. 

Fwake  ~  ^^CVmax^"max>  max^ lii/^^max)  (-.17) 

Uj,f  is  the  magnitude  of  the  velocity  profile's  velocity  range.  Fmtx  is  the  maximum  value 
provided  by  the  following. 

fW-j-MD  (2.18) 


is  the  value  of  y  corresponding  to  Fmtk.  Constants  in  the  Sankar  code  are  defined 
as  follows: 

A+~26.0 


k  —  0.4 


K  =  0.0168 

QP-  1-6 

FfCleb  =  0-3 

Cwk  =  0.25 

In  order  to  account  for  turbulence  outside  the  boundary  layer,  such  as  that  which 
occurs  in  the  dynamic  stall  process,  a  modified  turbulence  model  is  available  which  ef¬ 
fectively  increases  the  outer  model's  mixing  length.  In  this  case,  Fm„  and ymtl  are  deter¬ 
mined  by  the  following 

F(y)=>'2MD  (2.19) 


t 
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The  modified  turbulence  model  was  developed  in  response  to  early  predictions  of  How 
separation  at  high  angles  of  attack  and  has  been  found  to  cause  premature  reattachment 
during  the  downstroke  of  dynamic  cycles.  Its  use,  therefore,  is  only  recommended  until 
stall  onset  [Ref.  S.p.  33]. 

F.  CODE  STRUCTURE 

The  code  is  structured  such  that  the  following  user  defined  parameters  are  input 
from  logical  unit  05  (appended  to  the  end  of  the  code  for  Cray  processing):  grid  size,  step 
size  (dt),  artificial  viscosity  magnitude,  mean  angle  of  attack  (ot0),  oscillation  magnitude 
(a,),  angle  for  suspension  of  the  modified  turbulence  model,  reduced  frequency  (k).  free 
stream  Mach  number  (30,  Reynolds  number  (Re),  distance  of  the  first  n  ■  contour  from 
the  airfoil  wall,  starting  time,  pitch  and  restart  flags  and  airfoil  geometry  data.  Added 
to  this  list  for  the  present  study  are  the  number  of  time  steps  for  the  code  to  march  on 
the  present  run,  the  number  of  steps  between  each  plot  (output  interval)  and  the  total 
number  of  steps  and  plots  completed  on  all  previous  runs.  Variables  are  then  initialized, 
previous  stored  solution  datasets  are  retrieved  and  read  (to  allow  incorporation  of  all 
datasets,  including  those  from  the  present  run,  into  a  single,  combined  file)  and  flow  field 
starting  conditions  are  input.  A  series  of  subroutine  calls  then  generate  the  grid 
(AIRFOL,  SING,  TABINT,  WRAP),  cluster  grid  points  for  viscous  flows  (CLUSTR. 
STRTCH),  rotate  the  airfoil  (ROTGRID)  to  its  actual  angle  of  attack  and  compute  the 
initial  metrics  (METRIC). 

At  this  point,  the  code  is  fully  initialized  for  commencement  of  the  flow  field  sol¬ 
ution  process,  which  is  conducted  via  an  iterative  loop.  At  each  iteration,  time  is  first 
marched  forward  one  time  step,  followed  by  computation  of  the  time  dependent  values 
of  angular  velocity,  angle  of  attack  and  step  change  of  angle  of  attack,  according  to  the 
relations: 


co  =  2k\t  sin(2AA/^r) 


(2.20) 


a,-or0-ot,  co^lkM^i) 

«/_t  =  «o  “  «i  cosD/cd/^r  -  di)~] 
dot.  =  otj  -  oq_j 

The  grid  is  then  rotated,  in  the  physical  plane,  by  applying  the  the  following  relations 
at  each  grid  point. 


V'.'  -Z' 


/.a  A  A*\  AA -  lAA' 


x  =  x  cos  (c/a)  —  y  sin(c/a) 


(2-21) 


y  =y  cos(Ja)  -  x  sin(</sm'<?p/ta) 

Following  recomputation  of  the  metrics,  the  solution  is  computed  by  subroutine  SLPS, 
the  ADl-algorithm.  which  calls  DISSIP.  for  computation  of  the  explicit  dissipation 
terms,  STRESS  and  RESI,  for  computation  of  the  inviscid  and  viscous  terms,  respec¬ 
tively,  AMAT1  and  MATRIX1,  for  computation  of  the  Jacobian  and  its  inverse  in  the 
C  -direction  and  A.MAT2  and  MATR1X2,  for  computation  of  the  Jacobian  and  its  in¬ 
verse  in  the  rj  -direction.  Subroutine  STRESS  also  calls  subroutine  EDDY,  the  turbu¬ 
lence  model,  for  computation  of  the  viscosity  coefficient.  WALLBC  enforces  the  wall 
boundary  conditions  and  finally,  solution  files,  as  discribed  later,  are  generated. 

Prior  to  resumption  of  the  loop,  performance  coefficients  are  generated  by  subrou¬ 
tines  LOAD  and  CPPLOT.  Surface  pressure  coefficients  are  obtained  from  the  relation: 


Pb-P* 


WPoo  (-V4o«oo) 


(2.22) 


where  p„  is  the  surface  pressure  and  px  is  the  free  stream  pressure.  Skin  friction  coeffi¬ 
cients  are  a  function  of  the  wall  shear  stress  (r„)  according  to  the  relations: 


Cf- - — - -  (2.23) 

iv=uy-  vx  a  J{pnx;  +  y>{) 

The  aerodynamic  loads  of  lift,  drag  and  moment  about  the  quarter-chord,  are  computed 
by  the  following  relations. 


C,  =  C„  cos(a)  -  C,  sin(a)  (2.24) 

Q  =  Cn  sin(ot)  -  C,  cos(a) 

=  J  CpUy  -ye!A)dy  +  {x  -  xcji)dx] 
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C„  and  C,  are  the  normal  and  tangential  forces  obtained  from  summing  surface  pressure 
and  skin  friction  forces  according  to  die  following. 


C„  =  l  t>6r  + 


jc/v]/c 


(2.25) 


Cr  =  [  ~  J  C.Jy  +  J  C/.rJ/c 


G.  CODE  MODIFICATIONS 

Ihe  following  modifications  were  made  to  the  code,  resident  o^  the  FML  VAX. 
during  the  course  of  this  study. 

1.  Solution  output  commands  are  placed  within  the  loop  in  order  to  provide  interval 
output. 

2.  Read  commands  are  placed  at  the  beginning  of  the  main  program  in  order  to  allow 
storage  of  all  solutions  (previous  restarts  and  the  current  run)  in  a  single  combined 
file. 

3.  The  airfoil  section  is  rotated  to  the  initial  angle  of  attack,  vice  rotation  of  the  free- 
stream  direction. 

4.  Restarts  may  be  made  from  Plot3D  (output)  format  as  well  as  vector  format. 

5.  Additional  Plot3D  files  (surface  pressure  and  skin  friction  line  plots)  are  output 
from  subroutine  CPPLOT. 

6.  Grid  dimensions  are  made  true  variables. 

7.  U  ser  inputs  are  reformatted  for  ease  of  use. 
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III.  DYNAMIC  GRAPHICS  GENERATION 


A.  INTERACTIVE  COMPUTER  GRAPHICS 

Due  to  the  complexities  and  time  dependent  phenomenon  associated  with  unsteady 
Hows,  attainment  of  acceptable  clarity  in  flow  field  solutions  requires  complete  descrip¬ 
tive  information  across  fine  grids,  at  a  large  number  of  minute  time  intervals.  Thus, 
analysis  and  visualization  of  data  generated  by  codes  such  as  Sankar's,  is  difficult  due 
to  both  the  volume  of  information  provided  and  its  temporal  nature.  Requirements 
identified  for  effective  How  field  visualization  include  maximization  of  the  bandwidth  of 
information  transfer,  such  that  it  closely  matches  the  capabilities  of  the  human  eye, 
maximization  of  the  quality  of  graphical  displays,  by  enhancing  key  features  and  sup¬ 
pressing  others,  and  maximization  of  the  controllability  of  information  (Ref.  9|.  Addi¬ 
tional  advantages  in  information  transfer  can  be  achieved  through  the  use  of  redundant 
coding  and  structured  displays  [Ref.  10|.  In  response  to  these  requirements,  interactive 
computer  graphics  workstations  have  been  evolved  to  complement  the  super-computer. 
Workstation  capabilities,  in  terms  of  geometrical  transformation  and  screen  update  dis¬ 
patch  have  been  utilized  by  programmers  to  produce  effective  representations  of  flow 
field  motion,  through  synchronization  of  coordinated  data  sets.  Solution  clarity  is  en¬ 
hanced  by  the  high  degree  of  spatial  resolution  and  large  range  of  colors  afforded  by  the 
workstation  displays.  The  interactive  capabilities  which  currently  exist  for  semi-complex 
flows  serve  to  improve  visual  cues  for  display  of  three-dimensional  data  sets  and  allow 
immediate  access  to  regions  of  interest. 

B.  HARDWARE 

The  NASA-ARC  Fluid  Mechanics  Laboratory  is  currently  equipped  with  an 
IRIS-3000  series  workstation  configured  with  4M  bytes  of  display  list  memory  and 
equipped  with  z-clipping  and  z-buffering  hardware  and  a  multiple  key  mouse.  The  IRIS 
display  processing  unit's  refresh  memory  is  arranged  as  a  two-dimensional  array  (768  x 
1024),  with  row-column  positions  matching  display  pixel  x-y  coordinates.  24  bits  are 
reserved  for  each  pixel  at  eight  bits  (256  intensities)  per  color,  resulting  in  a  range  of 
16,777,216  possible  colors.  When  displaying  dynamic  graphics,  the  refresh  memory  is 
divided,  in  order  to  meet  user  demands,  since  flyback  time  (about  1.3  n  sec)  is  too  short 
for  complete  picture  update.  This  double  buffering  mode  utilizes  half  the  memory  for 
refreshing  the  screen  displays,  while  the  other  half  is  dynamically  updated.  This  ensures 


smooth  motion  on  the  display,  but  divides  the  number  of  bitplanes  per  pixel,  in  half. 
Thus,  the  number  of  colors  available  drops  to  4096.  In  order  to  avoid  a  similar  decrease 
in  the  range  of  colors  available,  color  maps  are  utilized.  In  this  case,  pixel  values  in  the 
refresh  memory'  are  routed  to  an  index  or  color  look-up  table  (with  each  entry  composed 
of  24  predefined  bits,  which  define  the  color),  vice  direct  routing  to  the  intensity  digital- 
to-analog  converter.  The  IRIS  transformation  rate,  via  a  single,  composed  transforma¬ 
tion  matrix  in  homogeneous  coordinates,  is  80,000  coordinates  per  second,  with  points 
and  lines  generated  at  3  million  pixels  per  second,  or  40.000  inches  per  second.  Polygons 
can  be  filled  (fiat  shading)  at  the  rate  of  44  million  pixels  per  second,  or  70,00<>  square 
inches  per  second.  These  capabilities  allow  generation  of  most  displays  at  the  rate  of 
approximately  10  screens  per  second,  which,  under  ordinary  conditions,  is  greater  than 
the  fusion  frequency  of  the  human  eye.  Interactive  transformations  are  ordinarily  ac¬ 
complished  via  the  mouse. 

C.  PLOT3D  SOFTWARE 

Solutions  provided  by  the  Sankar  code  consist  of  multiple  binary  files,  specifically 
formatted  for  input  into  the  Plot3D  graphics  software,  authored  by  Pieter  Buning  for 
NASA.  Two  filetypes  exist,  which  are  labeled  as  XYZ  or  grid  files  and  Q  or  flow- 
quantity  files.  Both  types  are  three-dimensional  matrices,  with  placement  of  data  within 
the  matrices  corresponding  to  the  row-column  addresses  of  mesh  points.  XYZ-files  have 
a  depth  of  two  in  the  two-'  imensional  case  and  provide  cartesian  coordinate  definitions 
of  each  mesh  point.  Q-files  have  a  depth  of  four  in  the  two-dimensional  case  and  pro¬ 
vide  computed  (low  quantities  at  each  mesh  point,  corresponding  to  the  q-vector  of 
equation  2.1.  Headers  for  each  file  list  free-stream  Mach,  Reynolds  number,  angle  of 
attack  and  time. 

Additional  flow  field  quantities  at  each  mesh  point  are  computed  according  to  the 
following  relations.  [Ref.  11] 

Definitions: 
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Based  on  these  values,  the  Plot3D  software  will  produce  plots  of  scalar  functions 
(density,  pressure,  temperature,  enthalpy,  energy,  velocity,  entropy,  momentum  and 
shocks)  suitable  for  contour  plots,  vector  functions  (velocity,  vorticity,  momentum  and 
pressure  gradient),  particle  trace  functions  (particle  traces  and  vortex  lines),  shock  wave 
locations  and  grid  functions  (computational  meshes  and  walls). 

Particle  traces  are  generated  using  trilinear  interpolation  of  values  of  the  vector 
function  inside  a  computational  cell  and  second-order  Runge-Kutta  steps  to  advance  the 
particle  in  space.  Five  steps  are  required  for  each  cell.  The  algorithm  for  computing 
shocks  computes  the  Mach  number  component  in  the  direction  of  the  local  pressure 
gradient.  Locations  where  this  value  decreases  through  1.0  is  plotted  as  a  shock.  Two- 
dimensional  stream  functions  are  calculated  by  integrating  the  mass  How  across  a  coor¬ 
dinate  line. 
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Output  from  Plot3D  is  available  in  a  variety  of  formats  and  is  a  function  of  user- 
defined  attributes.  For  dynamic  plots,  output  is  in  the  form  of  graphics  files,  formatted 
for  further  manipulation  by  animation  software.  Device  independent  plotting  (DIP)  is 
enabled  through  use  of  the  ARCGraph  binary  file  format  which  allows  data  to  be  ma¬ 
nipulated  bv  a  collection  of  function  libraries  and  utilities  developed  by  the  Advanced 
Computer  Research  Center  at  NASA-ARC.  Parameter  data,  along  with  graphics  prim¬ 
itives  (device  independent  op-codes)  arc  contained  in  these  files,  allowing  data  manipu¬ 
lation  by  the  Cray.  IRIS  and  attached  VAXes.  Plot3D  generation  of  graphics  files 
requires  that  several  attributes  be  defined.  In  order  to  reduce  user-tasking.  Plot3D  w  " 
initially  search  for  a  filename-specific  initialization  communications  file,  which  contains 
all  required  inputs.  Input  files  must  be  read  singularly,  which  makes  necessary  the  sep¬ 
aration  of  those  combined  plotting  files  received  from  the  Cray.  Output  graphics  tiles, 
however,  are  once  again  stored  in  combined  files  in  order  to  speed  further  processing. 
User-defined  attributes  include  axis  scales  and  extreme  values,  function,  number  of 
contours,  line  colors  and  types  and  wall  definitions. 

Subroutine  CPPLOT,  within  the  Sankar  code,  provides  grid  (XYZ)  files  which  con¬ 
tain  plotting  information  for  surface  pressure  coefficient  and  skin  friction  coefficient  line 
plots,  corresponding  to  each  output  flow  field  solution  file.  These  line  plots  are  scaled 
and  translated  after  separation  from  the  combined  files  and  plotted  utilizing  the  grid 
function.  Likewise,  an  angle  of  attack  pointer  plot,  utilizing  a  sine  wave  format  (gener¬ 
ated  external  to  the  code),  is  also  provided.  Q-files  containing  dummy  variables  (not 
utilized  for  grid  plot  generation)  are  also  provided  to  complement  the  line  plot 
XYZ-files,  as  required  by  Plot3D. 


D.  GRAPHICS  ANIMATION  SYSTEM  SOFTWARE 

The  final  step  in  the  visualization  process  is  the  animation  or  synchronization  of 
coordinated  graphics  datasets.  The  Graphics  Animation  System  (GAS)  is  a  software 
package  developed  by  Sterling  Software  under  contract  to  NASA  and  provided  to  edu¬ 
cational  institutions  throughout  the  country.  It  is  device  specific,  running  only  on  IRIS 
workstations,  and  is  written  in  the  C  progranuning  language  under  the  UNIX  operating 
system.  Graphics  files,  in  ARCGraph  format,  arc  read  and  stored  in  the  IRIS'  display 
list  memory  as  objects  by  the  menu-driven  program  and  assembled  according  to  se¬ 
quence  file  instructions  (stored  transformation  matrices).  In  order  to  smooth  motion 
associated  with  geometric  transformations,  up  to  30  linearly  interpolated  matrices  are 
provided  by  the  program,  between  user-identified  keyframes.  Titles,  legends  and  object 
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attributes  are  also  available.  Ultimately.  llo\v  field  representations  mi', 
to  video  tape  or  16mm  film.  Interactive  \iewiim  (oidinarilv  available  n 
controlled  via  a  mouse,  within  the  menu's  "view  data  "  window. 

A  complete  description  of  GAS  and  associated  videu  hardware  mu 
Reference  12. 
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IV.  RESULTS  AND  DISCUSSION 

A.  PROCEDURAL  CONSIDERATIONS 

Integration  of  the  FML  IRIS  with  related  CI:D  software  and  the  Cray  X-.V1P  48 
was  tested  by  plotting  the  flow  field  about  an  airfoil  experiencing  deep  dynamic  stall,  as 
provided  by  the  Sankar  code.  The  advantages  of  plotting  dynamic  stall  are  two-fold. 
First,  it  is  an  extremely  complex  and  time-dependent  phenomenon  which  is  difficult  to 
visualize  from  discrete  data  sets.  Thus,  the  accuracy  of  solutions  for  this  type  llow  field 
can  only  be  deternuned  by  consideration  of  the  complete  llow  field  in  both  time  and 
space.  Dynamic  graphics  are  therefore  ideally  suited  for  its  study.  Secondly,  the  FML 
is  currently  heavily  involved  in  studies  of  dynamic  stall,  which  include  verification  of  the 
Sankar  code.  Cray  output  is  thus  used  to  full  advantage. 

In  order  to  provide  smooth  animations,  an  extremely  large  number  of  plotting  sets, 
provided  at  equal  time  intervals,  must  be  generated  by  the  code  during  the  oscillatory 
cycle.  Since  access  and  transfer  of  this  data  is  required  several  times  during  the  graphics 
generation  process,  it  was  found  that  storage  of  datasets  in  combined  files  provided  the 
greatest  efficiency,  especially  in  transfers  between  the  Cray  and  IRIS.  Software  resident 
on  the  IRIS  is  then  employed  to  separate  the  combined  file  into  individual  datasets. 
(Embedded  in  these  programs  are  schemes  to  scale  data,  as  necessary  to  control  their 
ultimate  on-screen  positions,  allowing,  for  example,  placement  of  line  plots  in  a  specific 
corner  of  the  display  or  plotting,  for  comparison,  two  flow  fields  side-by-side.)  Individ¬ 
ual  datasets  are  automatically  assigned  names  which  correspond  to  those  contained  in 
Plot3D  initiation  files,  also  resident  on  the  IRIS.  These  files  currently  contain  com¬ 
mands  for  the  processing  of  50  datasets.  Thus,  the  combined  file  separation  and 
AR.Cgraph  file  production  (Plot3X)  process  must  be  completed  at  intervals.  The  use  of 
repetitive  file  names  during  this  process,  however,  increases  automation  to  such  an  ex¬ 
tent  that  only  a  limited  number  of  user  inputs  are  required  to  complete  the  process. 
ARCgraph  files  as  provided  by  Plot3D  are  again  in  combined  file  format,  allowing  easy 
(a  single  user  command)  input  into  animation  software. 

B.  THE  DYNAMIC  STALL  PROCESS 

A  prerequisite  to  review  of  any  code-provided  flow  field  solution  is  consideration  of 
available  empirical  information.  The  following  dynamic  stall  characteristics  have  been 
observed  for  both  harmonically  oscillating  airfoils  [Ref.  1 3 j  and  airfoils  undergoing 
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monotonic  angle  of  attack  increases  [Refs.  14,  15|.  First,  the  airfoil  stalls  at  an  angle  of 
attack  which  exceeds  the  static  stall  angle.  Second,  the  corresponding  normal  forces  and 
pitching  moments  exceed  those  attainable  in  the  static  case  and,  finally,  a  rapid  change 
in  pitching  moment  magnitude,  termed  as  "moment  stall"  by  Harris  and  Pruyn  [Ref.  16] 
occurs  several  degrees  in  azimuth  prior  to  the  normal  force  decrease,  or  "lift  stall",  vice 
simultaneous  occurence  as  in  the  static  case.  The  process  which  results  in  these  char¬ 
acteristics  can  be  broken  down  into  a  series  of  chronological  events,  as  depicted  in  Fig¬ 
ure  2  from  Reference  13.  While  the  specific  characteristics  for  a  given  airfoil  are  a 
function  of  free-stream  Mach  number.  Reynolds  number,  reduced  frequency,  mean  angle 
of  attack  and  oscillation  amplitude,  empirical  data  obtained  from  a  NACA-0012  airfoil 
at  k  =  .15,  Re  =  2.5xlO«  and  a  =  15’  —  10’  sin(a>/)  is  provided  in  the  following  discussion 
to  illustrate  the  dynamic  stall  process. 

At  point  (a)  of  Figure  2.  the  static  stall  angle  is  exceeded  without  any  detectable 
change  in  flow  over  the  airfoil.  The  boundary  layer  remains  thin  with  no  evidence  of 
flow  reversal  at  the  surface.  At  point  (b),  (  a  =  19  —  20”)  flow  reversal  appears  at  the 
surface.  Over  the  majority  of  the  airfoil,  the  boundary  layer  remains  thin  and  attached. 
At  the  rear  portion  of  the  airfoil,  however,  the  boundary  layer  thickens  as  the  flow 
gradually  decreases  to  zero  velocity.  At  point  (c),  large  eddies  appear  in  the  boundary 
layer.  By  point  (d),  flow  reversal  has  spread  over  much  of  the  chord,  from  the  trailing 
edge  to  x  c  =  0.3.  With  as  much  as  50%  of  the  airfoil  experiencing  flow  reversal,  no 
changes  in  C„  or  C„  are  detected.  At  point  (e),  (  a  =  23.4'  )  a  vortex  forms.  The 
boundary  layer  at  the  front  of  the  airfoil  abruptly  breaks  down  (simultaneously  from  x  c 
=  0.0  to  0.3)  and  the  vortex  moves  downstream  at  approximately  35  -  45%  of  the  free- 
stream  velocity.  At  point  (0  the  lift-curve  slope  increases  beyond  the  2rra  limit  of 
quasi-static  flows.  By  point  (g),  the  pressure  distribution  is  altered  sufficiently  to 
produce  a  noticeable  divergence  in  the  pitching  moment  (moment  stall)  but  is  accom¬ 
panied  by  a  continued  increase  in  lift.  Maximum  lift  occurs  at  approximately  mid-chord, 
followed  by  a  sharp  decrease  (point  (h)),  which  identifies  lift  stall.  (Boundary  layer  sep¬ 
aration  has  occurred  to  such  an  extent  that  continued  increases  in  angle  of  attack  will 
not  result  in  continued  increases  in  lift.)  At  point  (i)  (  a  =  24.95  °  ),  maximum  negative 
moment  occurs.  The  vortex  moves  off  the  trailing  edge  at  a  =  24.8  0  (downstroke)  and 
C„  and  C„  move  toward  static  levels.  At  point  (j),  the  airfoil  is  fully  stalled.  At  point  (,  k ), 
boundary  layer  flow  begins  to  reattach,  progressing  from  the  leading  edge  at  approxi¬ 
mately  25  -  35%  free-stream  velocity  and  is  complete  by  a.  =  7  0  (downstroke).  At  point 


(1),  after  some  hysteresis  (  a  =  6  0  ,  upstroke)  the  separated  region  has  fully  closed  and 
unstalled  values  of  C„  and  Cm  are  reestablished. 

C.  PROCEDURAL  VERIFICATION 

The  case  chosen  for  procedural  verification  closely  corresponds  to  the  Carr  data  of 
the  previous  section  and  exactly  matches  conditions  previously  run  when  the  code  was 
initially  vectorized  for  use  on  the  ARC  Cray  [Ref.  17).  There  are  several  advantages 
associated  with  this  duplication.  First,  it  allowed  storage  of  a  complete  record  of  this 
case  (data  saved  at  300  intervals,  vice  12)  for  future  access.  Second,  it  provided  a  simple 
means  by  which  to  determine  if  code  modifications  were  correctly  incorporated.  Finally, 
it  provided  a  baseline,  in  conjunction  with  empirical  data,  against  which  future  sensitiv¬ 
ity  studies  could  be  compared. 

Inputs  for  this  case  were  as  follows: 

Re  =  3.45x  106 
M„  =  .2S3 
dt  —  .005  sec. 

>1  -min  =  .00005 
co  =  .151 

Artificial  Viscosity  =  5 
Grid  size:  157  x  40 

Oscillatory  Cycle:  a  =  15°  -  10°  cos(co/) 

The  original  Baldwin- Lomax  turbulence  model  was  used  throughout  the  oscillatory 
cycle.  Figures  3  through  7  are  coefficient  of  pressure  contour  plots  of  the  predicted  flow 
field.  (Dynamic  plots  are  similar  to  these  Plot3D-provided  plots,  but  do  not  contain 
explicit  quantitative  information.)  As  previously  noted  [Ref.  17,  p.  48],  under  these 
conditions,  moment  coefficients  are  consistently  underpredicted  as  compared  to  exper¬ 
imental  data.  Drag  and  lift  coefficients  closely  match  experimental  data  below  approxi¬ 
mately  15  and  18  degrees  angle  of  attack,  respectively.  Early  prediction  of  flow 
separation,  however,  forces  inaccuracies  in  quantitative  solutions  beyond  these  values 
and  underprediction  of  the  maximum  lift  coefficient  obtainable.  These  findings  highlight 
the  rationale  behind  the  modified  turbulence  model.  Information  provided  by  the  dy¬ 
namic  plots,  shows  general  qualitative  agreement  with  experimental  data,  during  the 
upstroke.  A  vortex  forms  near  the  leading  edge  of  the  airfoil  section  and  progresses 
down  its  upper  surface,  gradually  increasing  in  size.  As  the  downstroke  begins,  however, 
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rather  than  shedding  from  the  trailing  edge,  the  vortex  stagnates  at  approximately  the 
65%  chord  point  and  gradually  dissipates.  Line  plots  of  surface  pressure  coefficients 
more  accurately  portray  (as  compared  to  traditional  carpet  plots)  the  turbulent  effects 
associated  with  the  existence  of  the  vortex.  During  the  vortex  dissipation  phase,  these 
line  plots  show  a  slight  (gradually  decreasing)  pressure  rise  in  the  vicinity  of  the  vortex. 
The  subtle  differences  between  these  plots  and  those  which  would  have  been  obtained 
had  the  vortex  actually  shed,  are  lost  when  this  information  is  integrated  lor  lift,  drag 
and,  to  a  lesser  extent,  moment  coefficient  data.  Thus,  quantitative  data  may  not  be 
good  indicators  of  flow  field  solution  accuracy. 

D.  SENSITIVITY  ANALYSIS 

The  following  code  parameters  which  should  effect  dissipation  of  the  vortex  have 
been  identified. 

1.  Artificial  viscosity  magnitude. 

2.  Turbulence  model  parameters. 

3.  Grid  resolution. 

A  limited  sensitivity  analysis  into  the  effects  of  grid  resolution  changes  in  the  trailing 
edge  region  was  conducted  in  order  to  provide  direction  for  follow-on  studies.  During 
the  initial  phase  of  this  analysis,  only  qualitative  results  will  be  considered,  for  those 
reasons  mentioned  above.  Dynamic  graphics,  therefore,  are  ideally  suited  for  this  task, 
especially  when  data  scaling  capabilities  are  utilized  in  order  to  plot  comparative  data¬ 
sets,  side  by  side. 

The  grid  is  relatively  coarse  in  the  region  in  which  the  vortex  becomes  stationary  and 
dissipates,  as  compared  to  those  regions  in  which  it  behaves  as  expected.  In  the  first 
sensitivity  case,  therefore,  grid  dimensions  were  held  constant  and  the  distance  of  the 
first  constant-  rj  line  from  the  airfoil  surface  was  increased  to  .001  ,  resulting  in  a  finer 
grid  in  the  trailing  edge  region.  Figure  S  shows  the  effects  of  this  change  on  the  vortex. 
Within  this  grid,  the  vortex  once  again  becomes  stationary,  but  effectively  splits,  such 
that  a  secondary  vortex  is  shed  while  the  original  vortex  dissipates.  While  this  does  not 
correspond  to  any  process  encountered  in  actual  flows,  it  does  indicate  that  final  sol¬ 
utions  are  sensitive  to  this  parameter. 

Changes  in  >i  -spacing  have  two  effects  on  the  flow  field  solution.  Grid  spacing  in 
the  trailing  edge  region  becomes  more  fine  at  the  expense  of  grid  spacing  in  the  boundary 
lacer,  which  becomes  more  coarse.  The  extent  of  influence  of  the  boundary  layer  sol- 


ution  on  the  vortex,  or  far  field  solution,  is  not  yet  clear.  It  is  known,  howeser.  that 
boundary  layer  solutions  are  extremely  sensitive  to  1/  -spacing  parameter  changes,  re¬ 
quiring  an  initial  value  of  .00005  for  accurate  solutions  (Ref.  IS].  Further  study,  there¬ 
fore.  is  required  in  order  to  determine  whether  vortex  behavior  changes  observed  in  the 
previous  case  were  in  response  to  the  altered  boundary  laser  solution  (which  would  re¬ 
commend  investigation  of  turbulence  model  parameters)  or  the  altered  grid.  In  order  to 
facilitate  this  study,  the  existing  grid  would  require  enlargement  (from  161  x  41).  Un¬ 
fortunately,  this  results  in  such  a  large  memory  allocation  from  the  Cray,  that  it  is  pro¬ 
hibitively  inefficient.  Replacement  of  the  entire  grid  generation  process  is  therefore 
currently  under  consideration  for  this  code. 

E.  CONCLUSIONS 

The  current  study  has  resulted  in  completion  of  the  following  items. 

1.  Full  procedural  documentation  (Appendix  A)  has  been  developed  for  efficient, 
code-independent,  two-dimensional,  dynamic  graphics  production  on  FML  and 
NFS  IRIS  workstations.  Data  must  be  in  PlotaD  format  and  may  be  provided  by 
either  code  or  experiment.  A  method  for  simultaneous  viewing  of  multiple  datasets 
has  been  incorporated. 

2.  The  Sankar  Navier-Stokes  solver  (Appendix  B)  has  been  modified  10  allow  dynamic 
graphical  presentation  of  its  flow  field  solutions. 

3.  Dynamic  graphics  applications  in  the  study  of  complex  flows  were  illustrated  by 
means  of  an  introductory  sensitivity  analysis,  using  the  Sankar  code.  This  limited 
analysis  identified  either  grid  resolution  in  the  trailing  edge  region  or  the  boundary 
layer  solution  as  parameters  to  which  vortex  behavior  is  strongly  dependent.  (This 
suggests  that  incorporation  of  a  more  sophisticated  grid  generation  scheme  would 
be  worthwhile.) 

In  order  to  derive  full  benefit  from  application  of  this  technology,  existing  and  future 
empirical  data  should  be  stored  in  plotting  format,  allowing  the  availability  of  visual  re¬ 
cords  of  multiple  flow  functions.  Also,  modification  of  additional  Navier-Stokes  solvers 
for  dynamic  output  will  allow  effective  code-to-code  and  code-to-e.xperiment  compar¬ 
isons,  all  in  an  effort  to  eventually  develop  an  accurate  flow  field  predictor. 


APPENDIX  A.  PROCEDURES  FOR  GENERATING  DYNAMIC 

GRAPHICS  AT  FML 


A.  THE  ANIMATION  PROCESS 

Creation  of  dynamic  graphics  tor  How  field  visualization  at  the  NASA-Amcs  Re¬ 
search  Center  Fluid  Dynamics  Laboratory  involves  five  basic  processes: 

1.  Navier-Stok.es  solver  is  submitted  to  the  Cray  X-MP  -IS  and  dynamic  output  saved 
in  the  Central  Storage  Facility  (CSF)  in  combined  files. 

2.  Combined  files  are  converted  to  IRIS  binary  and  transferred  to  an  FMLIRIS1  ac¬ 
count. 

3.  Combined  files  are  separated  into  individual  plotting  files  on  the  IRIS  and  con¬ 
verted  to  combined  graphics  (.gra)  files  by  Plot3D. 

4.  Graphics  files  are  fed  to  the  Graphics  Animation  System  software  (GAS)  and  dy¬ 
namic  graphics  are  plotted. 

5.  Resultant  files  and  plots  are  transferred  to  storage  or  video  tape. 

Prior  to  beginning  the  animation  process,  the  code  must  be  modified  to  provide  dy¬ 
namic  output  in  Plot3D  format.  This  involves: 

1.  Placement  of  solution  output  commands  within  the  iterative  loop  along  with  some 
interval  flag. 

2.  If  restarts  are  to  be  used,  placement  of  read  statements  at  the  beginning  of  the 
program  to  allow  storage  of  all  provided  solutions  in  a  single  combined  file. 

3.  Storage  of  any  additional  required  plotting  information  (line  plots)  in  Plot3D  XYZ 
files  for  eventual  display  utilizing  Plot3D's  grid  function. 

Full  documentation  of  Plot3D  formats  and  functions  may  be  found  in  Sterling  Software 
technical  note  15,  A  Guide  to  Generating  Movies  using  Plot3D  and  GAS,  by  Greg  Howe. 
The  procedures  outlined  below  refer  specifically  to  a  version  of  the  Sankar  Navier-Stokes 
solver,  modified  for  dynamic  output.  They  can  easily,  however,  be  generalized  to  apply 
to  any  code  (or  experiment)  from  which  dynamic  output  is  desired.  While  many  meth¬ 
ods  of  completing  these  steps  may  exist,  those  outlined  below  have  proven  to  be  the 
most  time  efficient. 

1.  Vax  Submitted  Cray  Jobs 

A  full  dvnamic  cycle,  utilizing  a  161  x  41  grid,  requires  over  5000  seconds  of 
Cray  time,  making  job  chaining  necessary  to  obtain  a  reasonable  priority.  With  900 
second  jobs,  24  to  4S  hours  are  usually  sufficient  for  the  full  computation  . '  be  com- 


pleted.  Since  significant  differences  exist  between  the  first  JCL  from  which  Plot3D  files 
are  saved  and  subsesquent  restart  JCL's,  two  JCL's  are  maintained  on  identical  copies 
of  Sankars  code. 

INITIAL. JCL.  Accesses  CSF  or  Cray  tape  (via  STAC j EX  commands)  for  initial 
solution,  which  may  be  in  either  Plot3D  or  matrix  format.  Combined  Plot3D  files  are 
initialized  on  CSF  and  assigned  PDVs  and  ID's.  Tape  <>5  input  data  (appended  to  end 
of  code)  is  assigned  in  accordance  with  definitions  included  in  MAIN  portion  of  the 
code.  (If  restart  is  from  dynamic  data.  AOA  will  not  be  a  whole  number.  TSTARf 
must  be  adjusted  accordingly,  in  order  to  have  accurate  time  AOA  matches.) 

RESTART.JCL.  Tape  05  inputs  must  identically  match  those  from 
IMTIAL.JCL.  with  the  exception  of  TSTART.CSTP.CPLT  and  possibly  FORMAT. 
INITIAL. JCL-dcfined  PDN  ID  solutions  are  accessed  (read,  stored  and  appended  with 
additional  soL.tions),  saved  (with  identical  PDN  ID's)  and  scrubbed  (oldest  versions 


deleted).  The  initial  restart  solution  is  accessed  with  librarv  routine  GLTP3D.  written 


by  Greg  Howe,  which  will  read  a  combined  file,  skip  a  specified  number  of  datasets 
(DSSKIP)  and  access  the  next  dataset  (or  datasets,  if  NL  MDS  is  not  equal  to  one,  the 
default).  Thus,  if 


DSSKIP  +  CPLT  -  1, 


the  correct  solution  will  be  provided  for  restarts.  Subsequent  restarts  illustrate  the  ad¬ 
vantage  of  combined  file  name  (PDN  ID)  repetition,  as  they  simply  require  the  following 
RESTART.JCL  and  Tape  05  updates: 


I.  DSSKIP 


2.  Job  chain  commands. 


3.  CS  TP 


4.  CPLT 


Once  the  optimum  number  of  plotting  sets  for  a  cycle  is  determined,  plotting 
files  for  an  AOA  pointer  on  IRIS  displays  can  be  generated  using  Plot3D.  (500  sets 
maximum,  or  modify  dimension  statements.) 

Line  Plots.  Any  pointer  or  line  plot  information  may  be  generated,  either  within 
the  code  or  externally,  and  stored  in  XYZ  files  for  further  processing.  Utilization  of 
Plot3D  function  0  and  proper  descriptors  of  walls  will  produce  the  desired  plots. 
SINE.TXT  is  a  program  which  generates  plots  for  a  dynamic  AOA  pointer  in  sine  wave 
format.  In  this  program.  Tape  05  includes  values  for  pointer  scaling  (axis  length  with 


mm 


Eg 


no  scalr  c  2 n  )  and  positioning  on  the  IRIS  screen  and  initial  pointer  luwuicn.  lYr 
examp.-.  '.he  cycle  starts  from  the  minimum  AOA. 

PNLOt.  I  =  ."(NPLOTS). 

Hardcopy  printouts  may  also  be  obtained  liom  these  grid  files  by  use  of  pic- 
grams  such  as  CPPI_OT.TXT.  This  program  provides  printouts  of  upper  and  lower 
surface  pressure  coefficient  and  skin  friction  values  at  \  c  locations,  including  tin  plot 
generation. 

2.  Path  to  the  IRIS 

SEND.JCL.  Utilizes  GETP3D  and  SENDWKS  to  convert  to  IRIS  binary  and 
transfer  complete  or  partial  combined  files,  or  individual  plotting  sets,  from  the  CRAY 
to  the  IRIS.  (Since  I  O  s  between  the  Cray  and  IR.  can  get  clogged'  if  too  many 
transfers  are  requested,  the  individual  method  is  not  recommended.)  This  may  need  to 
be  completed  in  stages  to  avoid  exceeding  the  IRIS'  memory  capacity. 

3.  Creating  Graphics  Files 

Combined  files  on  the  IRIS  are  separated  into  individual  plotting  sets  by  the 
programs  GETSIN.F  (also  generates  the  dummy  Q  file,  "qsin.iris ").  GETCP.F  (also 
scales  and  positions  Cp  or  skin  friction  on  the  IRIS  screen)  and  GETX.F  (for  XYZ  and 
Q  flow  field  datasets).  These  programs  require  some  initial  user  inputs  prior  to  running. 
Output  dataset  names  will  automatically  match  Plot3D  initialization  file  names.  Since 
the  PIot3x  run  for  each  dataset  type  (X.  P  and  SIN)  requires  specific  arguments,  initial¬ 
ization  files  for  each  data  type  exist  separately  in  Hies  XINT.COM,  PINT.COM  and 
SINT.COM.  Since  Plot3X  searches  for  the  file  name  PLOT3DINT.COM  ',  these  files 
must  be  renamed  appropriately,  utilizing  the  "mv"  .  unand.  Entering  tne  command 
"plot3x"  will  then  generate  multiple  graphics  files  wlu.n  are  stored  in  a  combined  file 
matching  the  initial  Q  file  name  encountered  by  Plot3D,  (i.e.,  "qOOl.gra").  This  file  is 
then  stored  on  an  appropriate  subdirectory  or  fed  to  GAS. 

4.  Notes  on  GAS 

For  movies,  the  maximum  number  of  objects  per  sequence  is  titty.  The 
"seq(n).seq"  files  are  50  object  far-field  solution  files  (seq(n).seq  =  objects  ni.louu, 
which  use  object  400  for  titles.  They  can  be  fed  to  GAS  (after  object  input)  via  the  AI  N 
1  O  window.  For  interactive  viewing,  the  VIEW  DATA  window  is  the  only  usable  op¬ 
tion.  (No  sequence  tiles  are  required  in  view  data.)  Full  GAS  documentation  is  avail¬ 
able  in  the  GAS  User's  Guide,  available  from  Staling  Software. 


5.  Storage 

Since  IRIS  memory  space  is  severely  limited  at  FML,  resultant  graphics  files 
should  be  stored  elsewhere,  if  not  in  immediate  use.  Downloading  to  1  4"  cassette  tape 


sily  accomplished  (also  a  recommended  backup),  as  is  transfer  to  CSf. 
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B.  SUPPORT  CODES 
1.  INITIAL.JCL 


JOS, JhNFLYNAVY.T-900.  ERIC  PAGENKOPF . *4269 . 

ACCOUNT , AC- . US- . UPW- . 

•  . 

CFT.ON-A.OFF-S. 

•  . 

•  RESTART  COMMANDS: - — - 

CSF . ACCESS . CSFACCT- . CSFPSWD- , DN-OLDSLN , PDNNNSE509 . ID-VALDES . 

ASSIGN . DN-OLDSLN , A-FT07 . 

•  .  —  “  —  —  "  '  ■  -  '  "  ■'  '  —  i  - 

..SAVE  DATA  FROM  PRESENT  RUN  IN  COMBINED  FILE: - 

ASSIGN.ON-xrZ.  A-FT20. 

ASSIGN. DN-Q.  A— FT21 . 

ASSIGN. DN-CPX ,  A-FT50 . 

ASSIGN, DN-CFX.  A-FT60. 

•  . 

LDR . MAP-PART , SET-ZERO . 

•  . 

CSF . SAVE , CSFACCT- . CSFPSWD- . 0»*XYZ .  PDN-SNKR05X , 

CS  F . SAVE . CSFACCT- . CSFPSWD- . DN-Q .  PDN-SNKR05Q . 

CSF .  SAVE .  CSFACCT- .  CSFPSWD- .  DN-CPX  .  PDN-SNKR05P  . 

CSF . SAVE . CSFACCT- , CSFPSWD- . DN-CFX .  PDN-SNKR05F , 

•  .  — —  I.  — — — — — — — — —  — 

• .ACCESS. DN-SENDVAX.PDN-SENDVAX. ID-STTRDM.OWN-RFTRDM. 

• .SENDVAX.DN=XYZ.VDN-’RAL"J1ANPS  poaaword" : : [ J I ANPS . PAGAN . DATA]XT1 .DAT’ . 

•  SENDVAX.DN-Q.VDN-,RAL"JIAh»>S  paaaword"  :  :[ J I ANPS . PAGAN. DATA) QT1  .DAT’  . 

•  ACCESS. DN-SENDWKS. PDN-SENDKKS. ID-STTRDM.OWKNRFTRDM. 

•  . SENDWKS . DN-Q3X49 . VDN- ' FMLIRIS1 “jiao  paasword" : : M/u/j i oo/pogon/q. i r i a*”  .NOREC. 

•  .  SENDWKS  ,  Df^IFX .  VDn-'  FMLIR1S1 "  j  ioa  poaaword'*:  :*'/u/j  i  oo/pagon/cf  01  dot'”  .NOREC. 

• .  ■—  -  — .  . . 

•JOB  CHAINING  COMMANDS: - 

FETCH . Dt^JOB2 . TEXT- ' [ PAGAN  JNSMULT . TXT ; 2 ' . 

REWIND. DN-JOB2. 

SUBMIT. ON-JOB2. 


ID-FLYNAVY. 

ID-FLYNAVY. 

ID-FLYNAVY. 

ID-FLYNAVY. 


2.  RESTART.JCL 


JOB. JN-FLYNAVY.T-900.  ERIC  PACENKOPF . X4269 . 

ACCOUNT  .  AO ,  US- .  UPW- . 

•  . 

CFT , ON— A , OFF— S . 

• . 

•RESTART  COMMANDS : - 

.  .  CSF . ACCESS . CSFACCT- . CSFPSWO . DN-OLDSLN , PON-NSE509 , 
CSF , ACCESS . CSFACCT- . CSFPSWO , DN-GETP3D .  PDN-GETP3D , 
CSF . ACCESS . CSFACCT- . CSFPSWO- , DN-XYZ 1 .  PON-SNKR05X , 

CSF .  ACCESS .  CSFACCT- .  CSFPSWO- .  DN-G 1  .  PDN-SNKR05Q . 

CSF, ACCESS, CSFACCT-. CSFPSWO-. DN-CPX1 ,  PDN-SNKR05P. 

CSF,  ACCESS. CSFACCT-. CSFPSWO. DN-CFX1  ,  PDOSNKR05F , 

REWIND. DN-Q1  . 

GETP3D . I -01 .OOLDSLN.DATATYPE-Q.DSSKIP-49. 

REWIND. DN-OLDSLN. 

REWIND. DNNXYZ1  . 

REWIND . DN-Q1 . 

REWIND . DN-CPX1  . 

REWIND.DN-CFX1  . 

ASSIGN. DN-XYZ 1 .  A-FT90 . 

ASSIGN . DN-Q1 ,  A-FT91. 

ASSIGN. DN-CPX1 .  A-FT92. 

ASSIGN . DN-CFX1 .  A-FT93 . 


IOVALDES. 

IORFRICC. 

I OFLYNAVY . 
IOFLYNAVY. 
2  OFLYNAVY. 
I  OFLYNAVY. 


ASSIGN. DN-XYZ 1 .  A-FT90 . 

ASSIGN . DN-Q1 ,  A-FT91. 

ASSIGN. DN-CPX1 .  A-FT92. 

ASSIGN . DN-CFX1 .  A-FT93 . 

ASS  I GN . DN-OLDSLN . A-FT07 . 

«  _  — —  '  -  ■  -n.i  —  -  i- 

•SAVE  DATA  FROM  PRESENT  RUN  IN  COMBINED  FILE- 
ASSIGN  .  DN-XYZ  .  A— FT 20 . 

ASSIGN, DN-0.  A— FT21 . 

ASSIGN. DN-CPX ,  A-FT50. 

ASSIGN.DN-CFX,  A-FT60. 


LDR.MAP-PART.SET- 

•  . 

CSF. SAVE. CSFACCT- 
CSF . SAVE .CSFACCT— 
CSF. SAVE. CSFACCT- 
C3F , SAVE . CSFACCT- 
CSF.  SCRUB.  CSFACCT. 
CSF,  SCRUB.  CSFACCT. 
CSF.  SCRUB.  CSFACCT. 
CSF.  SCRUB.  CSFACCT. 


.CSFPSWO.DANXYZ. 
.CSFPSWO  ,DN*Q, 
.CSFPSWO.  DN-CPX. 
.CSFPSWO. DNKIFX. 


PDN-SNKR05X, 

PDN-SNKR05Q. 

PDN-SNKR05P. 

PDN-SNKR05F, 


IOFLYNAVY. 
I  OFLYNAVY. 
IOFLYNAVY. 
IOFLYNAVY. 


. CSFPSWO. PDN-SNKR05F.  IOFLYNAVY. 
. CSFPSWO. PDOSNKR05P,  IOFLYNAVY. 
. CS FPSWO . PDOSNKR05X  .  IOFLYNAVY. 
. CSFPSWO. PDN-SNKR05O.  IOFLYNAVY. 


•  . ACCESS . DN-SENDVAX . PON— SENDVAX . I O .  OWN- . 

• . SENDVAX . DN-XYZ . VDN— • RAL" J I ANPS  password" :  :  {  J  I ANPS . PAGAN . DATA]XT 1 . DAT ’ 

•  . SEND VAX. DN-Q. VDN-’RAL"JJANPS  password": :[ J I ANPS. PAGAN. DATA JOT  1 . DAT ' . 

• . ACCESS . DN-SENDWKS . PDN— SENDWKS . I O . OWN- . 

• . SENDWKS .DN-Q3X49 .VDN-' FULIR I SI " j i aa  password" : : "/u/ j i oa/pogan/q . iris" 

•  . SENDWKS .DN-CFX .VDN-’ FMLIRIS1 " j i aa  password" : : "/u/j i ao/pogon/cf 01 .dot" 


. NOREC . 
,  NOREC . 


•JOB  CHAINING  COMMANDS: - 

FETCH . DN-J0B3 . TEXT- 1 ( PAGAN] 4063 . TXT ' . 
REWIND.  DOJ0B3 . 

SUBMIT . DN-J0B3 . 


,»4Vl 


.’♦’A*/!*  A* *•»*  1  Vlil’l 
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3.  SEND.JCL 


JOB, JN-GETX.T-30.  ERIC  PAGENKOPF,  X4269. 

ACCOUNT . AC- . US- . UPW- . 

•  . 

CSF . ACCESS . CSFACCT- . CSFPSWD-  . DN-XYZ . PDN-SNKR38X . I D-FLYNAVY . 

CSF , ACCESS . CSFACCT- , CSFPSWD- , DN-Q . PDN-SNKR38Q . I D-FLYNAVY . 

•  . 

CSF . ACCESS . CSFACCT- , CSFPSWD- , DN-GETP3D . PDN-GETP3D . I D-RFRICC . 

ACCESS . DN-SENDWKS . PDN-SENDWKS . I C*  .OWN- . 

•  . 

REWIND. DN-XYZ. 

REWIND. DN-Q. 

•  . 

GETP3D. I-XYZ , O-XOUT 1 , DATATYPE— XYZ .DSSK I P—1 00 , NUMDS— 1 50. 

GETP3D, I-Q.OQOUT1 , DATATYPE-0 , OSSK IP-1 00 , NUMDS-1 50 . 

•  . 

SENDWKS . DN— XOUT 1  ,  VON-  '  FMLIRIS1  "jiao  password"  :  :  *'/u/j  i  ao/pagan/38x .  i  r  i 
SENDWKS .DN-QOUT 1 , VON— ’ FMLIRIS1 “ j i oa  password" : : "/u/ j i aa/pogon/38q . i r i 

•  . 

• . FETCH. DN-J0B2. TEXT-1 [ PAGAN] SEND. JCL; 2 1 . 

• .REWIND. DN-J0B2. 

• . SUBM I T , DN- J0B2 . 


s"' .NOREC. 
s"’ .NOREC. 


non  o  ooo 


’JTY, 


LlUMiUUWlI 


-I.  SIi\E.TXT 


JOB, JN-FLYNAVY.T-30.  ERIC  PAGENKOPF . X4269 . 
ACCOUNT . AC- , US- , UPW- . 

•  . 

ACCESS  .  DN-SENDWKS ,  PD^SENDWKS  .ID-.  OWN- . 


CFT.0f*-A.0FF-S. 


ASSIGN .  DN-X . 
ASSIGN.ON—Q, 


A-FT50 . 
A-FT60 . 


LDR . MAP-PART , SET-ZERO . 

•  . 

SENDWKS.DN-Q.VDN-’FMLIRISI" j iao  PASSWORD": :/u/j i oo/pogan/qs i n . iris’ .NOREC. 
SENDWKS ,DN— X , VDN- ' FMLIRIS1 “ j i aa  PASSWORD”: :/u/j i ao/pagon/snk r05s ' .NOREC. 
CSF . SAVE . CSFACCT— . CSFPSWD- , DN-X , PDhNSNKR05S . I D-F  LYNAVY . 


PROGRAM  POINTER 


THIS  PROGRAM  GENERATES  PLOT3D  FILES  WHICH  WILL  PROOUCE  AN 
AOA  POINTER  WHEN  FUNCTION  0  IS  SELECTED.  XY2  FILES  ARE  SAVED 
ON  CSF  FOR  FUTURE  ACCESS  VIA  PROGRAM  GETSIN.JCL.  A  SINGLE  Q 
FILE  OF  PROPER  DIMENSION  IS  SENT  DIRECTLY  TO  THE  IRIS  WITH 
FILENAME  QSIN.IRIS. 

VARIABLES: 

NPLOTS :  TOTAL  NUMBER  OF  PLOTS  IN  ONE  CYCLE. 

PNLOCI :  INITIAL  INTEGER  LOCATION  OF  POINTER 
PNLOC :  INTEGER  LOCATION  OF  POINTER. 

X SCALE:  AXIS  SCALING  FACTOR 

Y SCALE:  SINE  WAVE  SCALING  FACTOR 

XPOSIT :  X-POSITIONAL  FACTOR  (SCREEN  LOCATION) 

YPOSIT:  Y-POSITIONAL  FACTOR  (SCREEN  LOCATION) 


INTEGER  PNLOC. GDIM. PNLOC I 

DIMENSION  SINX(3 . 500) ,SINY(3,500) .01 (3 . 500) .02(3 . 500) 
DIMENSION  03(3.500) .04(3,500) 

DATA  01 .02.03.04/6000.1 ./ 

...  INITIALIZE 
READ  (5,1) 

READ  (5.2)  NPLOTS. PNLOCI . XSCALE , YSCALE , XPOSIT, YPOSIT 
PI— AT AN ( 1 . )»4. 

GDIM  -  3 
PNLOC  -  PNLOCI 

DO  10  L-1. NPLOTS 

...  GENERATE  PLOT 3D  POINTER  FILE 

DO  20  N— 1 , NPLOTS+1 

XLOC  -  N«2. ‘PI/NPLOTS 
SINY(I.N)  -  (0.0). YSCALE  +  YPOSIT 
SINX(I.N)  -  (XLOC). XSCALE  +  XPOSIT 
S I NY ( 2 , N )  -  (SIN(XLOC))* YSCALE  +  YPOSIT 
SINX(2.N)  -  (XLOC). XSCALE  +  XPOSIT 
!0  CONTINUE 

PLOC  -  PNL0C.2. .PI/NPLOTS 
S I NY ( 3 , 1 )  -  (0.0)« YSCALE  +  YPOSIT 


I 

I 

$ 


w 


SINX(3 , 1 )  -  (PLOC).XSCALE  +  XPOSIT 
S I  NY ( 3 , 2 )  -  (SIN(PLOC) )»YSCAL£  +  YPOSIT 
S I  NX  (  3 . 2  )  -  (PLOC).XSCALE  +  XPOSIT 

#RITE(50)  GOIM.NPLOTS 

WRITE(50)((SINX(I, J).  1-1 ,3),J-1 , NPLOTS) , 

+  ((SINY(I.J).  1-1 .3) . J— 1 .NPLOTS) 

PNLOC  -  PNLOC  +  1 
I F( PNLOC . EQ . NPLOTS) THEN 

PNLOC  -  PNLOC  +  1  -  NPLOTS 

END  IF 

0  CONTINUE 

...  GENERATE  DUMMY  0  FILE 

WRIT£(60)  GOIM.NPLOTS 

WRITE(60)  GDIM.GDIM.GDIM.GDIM 

WRITE(60)  ((01(1. J).  1-1 .3) ,J-1 .NPLOTS). 

+  ((02(1. J).  1-1 .3). J-1 .NPLOTS). 

+  ((03(1. J).  1-1 .3), J-1 .NPLOTS). 

+  ((04(1. J).  1-1 .3) .J-1 .NPLOTS) 

FORMAT  (IX) 

FORMAT  (1X.2I8.4F8.4) 

END 

EOF 

NPLOTS:  PNLOC I :  XSCALE:  YSCALE :  XPOSIT:  YPOSIT: 

294  222  .08  .1  -1.0  -.5 


ms* 


5.  CPPLOT.TXT 


JOB . JN-SKYHAWK , T-30 . 
ACCOUNT .  AC- .  US- .  UPW- . 
*  . 

CFT.CXW-A.OFF-S. 


ERIC  PAGENKOPF , X4269 . 


CSF . ACCESS , CSFACCT- . CSFPSWD- , DN-GETP3D , PDN-CETP3D , I D-RFR ICC . 
ACCESS . DN-SENDWKS , PDMSENDWKS . I D-STTRDM . OWN— RFTRDM . 


PDN— SNKR05P , ID— FLYNAVY . 
PDMSNKR05F.  ID— FLYNAVY . 
PDf*SNKR05Q . I D-FLYNAVY . 


CSF. ACCESS, CSFACCT-. CSFPSWD-. DNW3PX1 .  PDN— SNKR05P 
CSF. ACCESS. CSFACCT-. CSFPSWD- , DN-CFX1 .  PDMSNKR05F 
CSF. ACCESS. CSFACCT-, CSFPSWD-. DN-Q 1 .  PD^SNKR05Q 

REWI NO . DN-CPX1 . 

REWI ND . DN-CFX1 . 

REWIND. DN-01  . 

GETP3D . I -CPX 1 . O-CPX . DATATYPE-XYZ . DSSK I  P-49 , NUMDS-2 . 
GETP3D . I-CFX1 ,  O-CFX, DAT AT YPE-XYZ, DSSK I  P-49 .NUMDS-2 . 
GETP3D . I-OI ,  CM3.  DATATYPE-Q.  DSSK IP-49 , NUMOS-2 ■ 
REWIND . DN-CPX . 

REWI ND  .  DN-CFX . 

REWIND. DN-Q. 

ASS  I GN . DN-CPX . A-FT50 . 

ASS I GN . DN-CFX . A-FT60 . 

ASSIGN. DN-Q.  A-FT70. 

•  . 

LDR . MAP-PART . SET-ZERO . 


PROGRAM  CPP LOT 


THIS  PROGRAM  READS  STORED.  TRUE  VALUES  OF  XX.  CP  AND  CF. 
AND  PRINTS  THEM  OUT  TO  TAPE  06. 

VAR 1 ABLES * 

NUMDS:  NUMBER  OF  OATA  SETS  IN  THE  READ  FILE 
NPLOT :  PLOT  NO.  OF  FIRST  DATA  SET  (DSSKIP  +  1) 


TAPE  ASSIGNMENTS: 


TAPE  05 
TAPE  06 
TAPE  50 
TAPE  60 
TAPE  70 


INPUT  DATA 
OUTPUT  DATA 
TRUE  CP  VALUES 
TRUE  CF  VALUES 
Q  FILES 


DIMENSION  CPX (3. 49) ,CFX(3,49) ,CPY( 3 , 49) ,CFY(3 . 49) 
DIMENSION  01(161 .41 ).Q2( 161 .41 ),Q3( 161, 41), Q4(161 .41) 
DIMENSION  K0DE(4) . LINE(90) 

DATA  K00E/1H  . 1H+, 1HI , 1H«/ 

READ  INPUT  DATA  k  INITIALIZE 

READ  (5.1) 

READ  (5.2)  NUMDS. NPLOT 
DO  9  1-1 .90 

LINE(I)  -  KOOE( 1 ) 

CONTINUE 

READ  k  SAVE  TRUE  VALUES 


O  *  Oi  M 


DO  10  LI-1 . NUMDS 

READ(50)  IMAX.KMAX 

READ(50)  ((CPX(I.J).  1-1 . IMAX) . J-1 .KMAX). 

+  ((CPY(I.J).  1-1 , IMAX) .J-1 , KMAX) 

READ (60)  IMAX.KMAX 

READ (60)  ((CFX(I.J).  1-1 . IMAX) . J-1 , KMAX) . 

+  ((CFY(I.J).  1-1 . IMAX) .J-1 , KMAX) 

READ (70)  IMAXO.KMAXO 

READ (70)  AMINF.ALFAD.REYREF.TIME 

READ (70)  ( (Q1 ( 1 , J ) ,  1-1 . IMAXQ) , J-1 .KMAXO) . 

+  ( (02 ( I . J ) ,  1-1. IMAXQ), J-1, KMAXO). 

+  ((03(1. J).  1-1 . IMAXQ) .J-1 .KMAXO) , 

+  ((04(1. J).  1-1 . 1MAX0) , J-1 .KMAXO) 

CP0  «  (1.  +  .2  •  AMINF..2)  ••  3.5  -  1. 

CP0  -  CP0  /  (  .7  •  AM1NF»»2) 

K0  -  30.  •  CP0  +4.5 
00  11  L-1.KMAX 

K  -  30.  •  (CP0  -  CPY{ 3 . L) )  +4.5 
K1  -  30.  •  (CP0  -CPY(2 ,  L)  )  +  4.5 
IF(K.LT.I)  K  -  1 
I F(K1 . LT . 1 )  K 1  -  1 
I F(K . GT . 90)  K  -  90 
I F ( K 1  .GT.  90)  K1  -  90 
LINE(K0)  -  K0DE(3) 

LINE(K)  -  K00E(2) 

LINE(KI)  -  K0DE(4) 

I F ( L . £0 . 1 )THEN 

WRITE(6.*)’  PLOT  AOA  TIME  MACH  REY  NO.  * 

WRITE(6.3)  NPLOT . ALFAD.TIME.AMINF.REYREF 
WR ITE(6 . • ) '  XOC  CFL  CFU  CPL  CPU  ' 

WRI TE(6 , 4)  CPX(1 .L).CFY(2.L).CFY(3.L).CPY(2,L).CPY(3.L).LINE 
ELSE 

V*RITE(6.4)  CPX(1  .L).CFY(2.L).CFY(3.L).CPY(2.L).CPY(3.L).LINE 
END  IF 

LINE(K1 )— KODE( 1 ) 

LINE(K)  -KOOE(I) 

11  CONTINUE 

NPLOT  -  NPLOT  +  1 
WRI TE(6 , 1 ) 

10  CONTINUE 

C 

1  FORMAT  (IX) 

FORMAT  (IX. 217) 

FORMAT  (IX. 14.3F9.4.F10.0) 

FORMAT  (1X.F6.4.4F8.4.90A1) 
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6.  GETX.F 


PROGRAM  GETX 
C 


THIS  PROGRAM  READS  COMBINED  (MULTIPLE  DATA  SET)  FILES 
OF  XYZ  AND  0  PLOTTING  DATA  AND  SEPARATES  THEM  INTO 
INDIVIDUAL  FILES  FOR  PLOT3D  SUBMITTAL. 

VARIABLES: 

CFXNAME:  NAME  OF  COMBINED  XYZ  FILE 
CFONAME :  NAME  OF  COMBINED  0  FILE 
FXNAME:  NAME  OF  INDIVIDUAL  X  FILES.  CHAR  VARIABLE 
FQNAME:  NAME  OF  INDIVIDUAL  0  FILES.  CHAR  VARIABLE 
DSSKIP:  NUMBER  OF  DATA  SETS  TO  SKIP  WHEN 
READING  COMBINED  FILE 
DSREAD :  NUMBER  OF  DATA  SETS  TO  READ 
INTVL :  INTERVAL  BETWEEN  DATA  SETS  SENT  TO  FILE 


C 

INTEGER  DSSKIP. DSREAD. COUNT 

CHARACTER  FXNAME* 1 2 . FQNAME • 1 2 .CFXNAME* 10 . CFQNAME* 1 0 
DIMENSION  X(250. 100) ,Z(250 . 100) 

DIMENSION  Q1 (250. 100). 02(250.100). 03(250, 100) .04(250.100) 
C 

•  ••  USER  INPUT:  ••• . .....«» . ••••• . . . 

C 

DATA  XSCALE . YSCALE/1 .  .1  ./ 

DATA  XPOS I T , YPOS I T / — 1 . . - 1 . / 

DATA  DSSKIP/10/.OSREAD/3/. INTVL/1/ 

CFXNAME  -  ' 38x .iris' 

CFQNAME  -  ' 38q .iris' 

C 


C 


C 


C 


10 

C 


OPEN ( UN I T-9 1 , F I LE-CFXNAME . STATUS- ' OLD 1 , FORM- ' B I  NARY • ) 
OPEN  ( UN  I T-92  .  F I  LE*=C FQNAME .  STATUS- '  OLD ’.  FORM- ’  B I  NARY  ’  ) 

FXNAME  -  '*000. iris' 

FONAME  -  ' q000 .iris' 

COUNT  -  0 


IF  (DSSKIP  GT.0)  THEN 

DO  10  ISKIP  -  1 .DSSKIP 


+ 

+ 

+ 


CONTINUE 


READ(91)  I  MAX , KMAX 

READ(9 1 )  ((X(I .K) . 1-1 . IMAX) ,K-1 ,KMAX ) . 

((Z(I.K).I-I . IMAX ) .K— 1 , KMAX ) 
READ (92)  IMAX. KMAX 
READ(92)  A.B.C.D 

READ(92)  ((01(1 .K), 1-1 . IMAX ) ,K— 1 .KMAX) , 
((02(1 .K) . 1—1 .IMAX) . K— 1 .KMAX) , 
( (03( I .K) . 1-1 . IMAX) ,K-1 .KMAX) , 
((04(1  K) . 1-1 , IMAX) ,K-1 .KMAX) 


END  IF 


DO  20  I  READ  -  1 .DSREAD 
COUNT  -  COUNT  -f  1 
WRITE  (FXNAME(2  4). 100)  IREAD 
WRITE  ( FONAME ( 2  4 ) , 1 00 )  IREAD 
READ(91 )  IMAX, KMAX 

READ( 9 1 )  ((X( I .K). 1-1 . IMAX) ,K— 1 .KMAX). 
+  ( (Z( I ,K) , 1-1 . IMAX) .K-l .KMAX) 


n 


•J 


3 
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12 

11 


20 

C 

100 


If  (COUNT. EQ. INTVl)  THEN 

OP£N(UNIT«1 , FI lE«FXNAME.FOR»M'B INARY ‘ ) 
DO  11  1-1 . I WAX 

DO  12  K-1 .KMAX 

X(  1  . K )— X ( I . K ) *X$CALE+XPOS I T 
Z(i  .K)-Z(I ,K).YSCALE+YP0S1T 
CONTINUE 
CONTINUE 

WRITE(I)  I MAX .KMAX 

WR1TE(1)  ((X(I ,K) . 1-1 . IMAX) . K-1 ,KUAX). 

( (Z( I .K) . 1-1 . IMAX) . K-1 .KMAX) 

CLOSE(I) 

END  IF 

READ(92)  IMAX. KMAX 
READ( 92)  A.B.C.D 

READ(92)  ((01(1 .K). 1-1 . IMAX). K-1, KMAX). 

((02(1 ,K) . 1-1 . IMAX) .K-1 .KMAX) . 

( (Q3( I ,K) . 1-1 , IMAX) ,K-1 .KMAX) . 

((Q4( 1 ,K) . 1-1 . IMAX) .K-1 .KMAX) 

IF  (COUNT. EQ. INTVl)  THEN 

OPEN ( UN I T-2 . F 1 L  E-FONAME , FORM*  -BINARY-) 
WRITE(2)  IMAX. KMAX 
*RITE(2)  A.B.C.D 

WRITE(2)  ( (01 ( I, K). 1-1 .IMAX). K-1 .KMAX). 

((02(1 .K). 1-1 .IMAX) . K-1 .KMAX). 
((03(1 .K). 1-1 .IMAX) , K-1 .KMAX). 
( (04(1 .K) , J-1 . IMAX) ,K-1 .KMAX) 

CL0SE(2) 

COUNT  -  0 

END  IF 
CONTINUE 


CL0SE(91 ) 

CL0SE(92) 

FORMAT (13.3) 

STOP 

END 
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7.  GETCP.F 


PROGRAM  GETCP 
C 


THIS  PROGRAM  ACCESSES  COACINED  CP  OR  CF  FILES  ON  THE 
IRIS  ACCOUNT.  SEPARATES  THEM  INTO  INDIVIDUAL  PLOTTING 
FILES  AND  SCALES  THEM  TO  MEET  SPECIFIED  REQUIREMENTS. 

A  DUfcMY  0  FI Lr  IS  ALSO  GENERATED  (Q.IRIS). 

VARIABLES: 

CPNAME:  COMBINED  FILE  NAME 

FNAME:  OUTPUT  FILE  NAME.  CHARACTER  VARIABLE 
DSSKIP:  NUMBER  OF  DATASETS  IN  THE  COMBINED  FILE 
TO  SKIP 

DSFILE:  NUMBER  OF  DATASETS  IN  THE  COMBINED  FILE 
TO  SEPARATE  AND  FILE 
X SCALE.  AXIS  SCALING  FACTOR 

YSCALE :  CP/CF  SCALING  FACTOR  (RELATIVE  MAGNITUDE) 
XPOSIT :  X  POSITION  FACTOR  (SCREEN  LOCATION) 
YPOSIT:  Y  POSITION  FACTOR  (SCREEN  LOCATION) 


C 

INTEGER  DSSKIP. DSFILE 
CHARACTER  FNAME* 1 t . QNAME.6 , CFNAME* 10 
DIMENSION  X(3.75).Z(3.75).XS(3.75),ZS(3,75) 
DIMENSION  01(3.75) .02(3 .75) .03(3 . 75) .04(3 . 75) 
DATA  01 ,02.03.04/900*1 ./ 

C 

•  ••  USER  INPUTS:  . . ••• . . 

C 

DATA  DSSKIP/250/, DSFI LE/44/ 

DATA  XSC ALE/. 5/. YSCALE/-. 05/ 

DATA  XP0SIT/-1 .0/. YPOSIT/. 5/ 

CFNAME  -  ’  ante  r05p  ’ 

FNAME-  'cp000. iris' 

C 


C 

OPEN( UNI  T-90.  FI  LE-CFNAME . STATUS- '  OLD FORInM '  BINARY ’  ) 

C 

QNAME  -  ' q . i r i s ' 

C 

«»•  SKIP  DSSKIP  FILES 
C 

IF  (DSSKIP. GT.0)  THEN 

DO  10  ISKIP  -  1 .DSSKIP 
READ (90)  I MAX . KMAX 

READ (90)  ( (X( I ,K) , 1-1 , IMAX) .K-1 .KMAX) , 

+  (<Z(I ,K) , 1-1 . IMAX) , K-1 .KMAX) 

10  CONTINUE 

END  IF 
C 

•••  SEPARATE  AND  SCALE  DSFILE  FILES 
C 

DO  40  IFILE  -  1. DSFILE 

WRI TE( FNAME(3 : 5) ,100)  IFILE 


» 

“A 


OPEN( UNIT-1  ,  FI  LE-FNAME.FORI»M' BINARY '  ) 
READ (90)  IMAX. KMAX 

READ(90)  ( (X( I .K) . 1—1 . IMAX) , K-1 .KMAX) . 
+  ((Z(I .K) . 1-1 . IMAX) , K-1 .KMAX) 
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DO  30  I  -  1 , I MAX 

DO  20  K  -  1 , KMAX 

XS( I , K )— X ( I ,K).XSCALE  +  XPOSIT 
ZS(I.K)— Z(I,K) • YSCALE  +  YPOSIT 
CONTINUE 
CONTINUE 

WRITE(I)  I MAX . KMAX 

#RITE( 1)  ((XS(l.K). 1-1 . IMAX) ,K-1 .KMAX). 

+  ((ZS( I ,K) . 1-1 . IMAX) .K-1 ,KMAX) 

CLOSE( 1 ) 

CONTINUE 

GENERATE  DUM^Y  0  FILE 

OPEN  ( UN  I  T-1  .  F I  LE-ONAME  ,  FORM-  '  8 1  NARY  '  ) 

WRITE(I)  IMAX. KMAX 

WRITE(l)  IMAX. IMAX. IMAX, IMAX 

WRI TE( 1 )  ((01 ( I ,K) , 1-1 . IMAX) , K-1 .KMAX) . 

+  ((02(1 ,K) , 1-1 , IMAX). K-1 .KMAX). 

+  ((03(1 ,K). 1-1 . IMAX), K-1 , KMAX ) . 

+  ((04(1 .K) , 1-1 , IMAX) .K-1 ,KMAX) 

CLOSE(I) 

FORMAT( 13.3) 

STOP 

END 


8.  GETSIN.F 


PROGRAM  GETSIN 


THIS  PROGRAM  READS  A  COMBINED  FILE  OF  AOA  POINTER 
DATA  (GENERATED  WITH  PLOT3D  FUNCTION  0)  AND  SEPARATES 
IT  INTO  INDIVIDUAL  PLOTTING  FILES. 


VARIABLES: 
CFNAME 
FNAME : 
DSSKIP 

DSFILE 


NAME  OF  COMBINED  FILE 
NAME  OF  INDIVIDUAL  FILES.  CHAR  VARIABLE 
NUMBER  OF  DATA  SETS  TO  SKIP  WHEN 
READING  COMBINED  FILE 

NUMBER  OF  DATA  SETS  TO  SEND  TO  INDIVIDUAL 


INTEGER  DSSKIP, DSFILE 
CHARACTER  FNAME* 1 2 .CFNAME* 1 0 
DIMENSION  X(3 , 294) , Z(3 , 294) 

USER  INPUT:  . . . 

DATA  DSSK I P/250/ .DSFILE/ 44/ 
CFNAME  -  ’ankrOSs- 


OPEN(UNIT-90. FI LE-CFNAME. STATUS- ’OLD’ , FORU- ' B I NARY  1 ) 

FNAME  -  • 9 i n000 . iris' 

IF  (DSSKIP. GT.0)  THEN 

DO  10  1SKIP  -  1 .DSSKIP 
READ(90)  I MAX , KMAX 

READ( 90)  ((X(I ,J). 1-1 . I MAX ) , J— 1 .KMAX). 

((Z(I.J).I— 1.1 MAX ) , J— 1 , KMAX ) 

CONTINUE 


DO  20  IFILE  -  1 .DSFILE 

WRITE  (FNAME(4.6) ,  100)  IFILE 

OPEN ( UN I T- 1 . F I LE-FNAME , FORM- ’ B I NARY ' ) 

READ(90)  I MAX . KMAX 

READ(90)  ((X(I.J).I  — 1,1  MAX) , J— 1  , KMAX) . 

((Z(I .J) , 1-1 . I MAX) , J-1 .KMAX) 
WRITE(I)  I MAX, KMAX 

WRITE(I)  ((X(I.J).I-I . IMAX).J-1 .KMAX), 

( (Z( I , J ) . 1  —  1 . I  MAX) . J-1 , KMAX) 

CLOSE(I) 

CONTINUE 

CLOSE(90) 

FORMAT (13.3) 

STOP 

END 
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C.  PROCEDURAL  DOCUMENTATION 


After  proper  initialization  (user-inputs)  as  described  abuse,  all  codes  or  J(  I.  >.  .'.iiu.ii 
produce  or  transfer  PlotaD  files  ate  submitted  to  the  (  ray  utilizing  stundai  J  (  SI  B 
commands.  Once  the  data  is  resident  on  the  IRIS  and  the  user  is  lotted  on  and  in  the 
proper  directory  level,  the  followin'.:  procedures  may  be  utilized  lor  graphics  generation. 

1.  Separate  the  combined  file  into  discrete  data  sets. 

Edit  "user  inputs"  section  of  GE  IX1'  (or  other  liics  as  appropriate)  to  d  ame  scal¬ 
ing  (ultimate  location  on  the  display),  and  identity  the  combined  lile  to  access  and 
those  data  sets  to  be  separated: 

vi  getx.f 

Compile  the  editted  program: 

177  getx.f 

Run  the  compiled  version  by  entering  the  filename  (a. out  by  dclauit). 

a. out 

Discrete  data  sets  will  appear  on  the  account  with  names  corresponding  to  those 
in  the  Plot3D  initialization  files.  Check  by  using  the  Is  command. 

2.  Use  Plot3D  to  generate  .gra  files. 

Use  the  mv  commands  to  give  the  proper  initialization  file  the  filename 
"plot3dini.com".  (The  IRIS  will  not  save  multiple  versions  of  files  with  the  same 
filename.  Instead,  older  versions  of  the  file  will  be  deleted.  In  ordo"  'o  u-.oid  in¬ 
advertent  deletion  of  initialization  files,  multiple  mv  commands  may  cc  required.) 

mv  xini.com  plot3dini.com 

If  processing  flow  field  files  (X  and  Qi.  edit  the  initialization  lile  to  identify  desired 
function  and  number  of  contours.  This  involves  changing  only  two  lines  at  the  top 
of  the  file. 


vi  plot3dini.com 

Create  .gra  files. 

plot3x 

Rename  the  resultant  combined  graphics  file. 

mv  qOOl.gra  somefilnm.gra 

After  creation  of  as  many  .gra  files  as  required  from  the  separated  liics.  clean  up 
the  account  using  the  following  wildcards. 


3.  Interactive  viewing. 

Run  the  Graphics  Animation  System  software. 

gas 

Input  .gra  files  to  GAS  by  selecting  on  the  main  menu: 

ARCGRAPH  file  input 

On  the  submenu,  select: 

load  entire  lile 

In  response  to  prompts,  enter  the  seoueutial  object  number  ami  !i,e  name 
( '  someilnin.gra  b.  For  color  plots,  press  UPTURN  when  prompted  for  tin: 
colormap. 

Mew  the  data  by  selecting  on  the  main  menu: 

view  data 

Utilize  the  mouse  to  manipulate  the  display. 


SANKAR  NAVI ER  STOKES  CODE:  SOLVES  TWO  DIMENSIONAL  FLOWS  PAST 
ARBITRARY  GEOMETRIES  USING  ADI  PROCEDURE.  THIS  VERSION  SAVES 
MULTIPLE  PLOT3D  DATA  SETS  IN  SINGLE  FILES. 

VARIABLES: 

TITLE:  (60  CHARACTERS  MAX. ) 

I  MAX:  NUMBER  OF  X  COORD  LOCATIONS 
KMAX:  NUMBER  OF  Y  COORD  LOCATIONS 
ITEL:  GRID  POINT  AT  LOWER  TRAILING  EDGE 
ITEU :  GRID  POINT  AT  UPPER  TRAILING  EDGE 
DT:  SIZE  OF  TIME  STEP 
WW:  EXPLICIT  ARTIFICIAL  VISCOSITY  TERM 
ALFA:  MEAN  ANGLE  OF  ATTACK  (DEGREES) 

ALFA1 :  AMPLITUDE  OF  OSCILLATION  (DEGREES) 

ALFAI:  ANGLE  BELOW  WHICH  MODIFIED  TURBULENCE  MODEL  USED  TO 
COMPUTE  EDDY  VISCOSITY 
REDFRE:  REDUCED  FREQUENCY 
AMINF:  FREE  STREAM  MACH  NUMBER 

REYREF:  REYNOLDS  NU»«ER  (MILLIONS;  NEG.  «  INVISCID  FLOW) 

DNMIN:  DISTANCE  OF  FIRST  POINT  OFF  OF  WALL 

TSTART :  TIME  FOR  CALCULATIONS  TO  START  ON  PRESENT  RUN 

FORMAT:  OLDSLN  FORMAT ; 3 . 0-PLOT 3D  FILES. -3. 0-0  MATRICES 

RSTRT :  TRUE  -  STORED  DATA  READ  TO  CONTINUE  ITERATION 

PITCH.  TRUE  -  AIRFOL  AOA  OSCILLATES 

PLUNGE:  TRUE  -  AIRFOIL  OSCILLATES  VERTICALLY 

FNU :  NUMBER  OF  UPPER  SURFACE  DEFINITION  POINTS 

FNL :  NUMBER  OF  LOWER  SURFACE  DEFINITION  POINTS 

FSYM:  SYMMETRY  FLAG  (1  -  SYMMETRIC) 

X:  AIRFOIL  GEOMETRY  DEFINITION  POINTS 
Y:  AIRFOIL  GEOMETRY  DEFINITION  POINTS 
CSTP:  NUMBER  OF  STEPS  COMPLETED 
CPLT:  NUMBER  OF  DYNAMIC  PLOTS  COMPLETED 

NSTP :  NUMBER  OF  TIME  STEPS  TO  BE  COMPLETED  ON  PRESENT  RUN 
PSTP:  NUMBER  OF  TIME  STEPS  BETWEEN  PLOT  DATA  OUTPUT 


TAPE  DEFINITIONS: 

TAPE  05:  INPUT  DATA 
TAPE  06:  OUTPUT  DATA 

TAPE  07:  FLOW  FIELD  INPUT  DATA  FOR  RESTARTS 
TAPE  20:  PLOT 3D  XYZ  FILE  STORAGE 
TAPE  21:  PLOT 3D  0  FILE  STORAGE 
TAPE  50:  PL0T3D  CP  XYZ  FILE  STORAGE 
TAPE  60:  PL0T3D  CF  XYZ  FILE  STORAGE 
TAPE  90:  PREVIOUS  RUN  X  FILE  STORAGE 
TAPE  91:  PREVIOUS  RUN  Q  FILE  STORAGE 
TAPE  92:  PREVIOUS  RUN  CP  FILE  STORAGE 
TAPE  93:  PREVIOUS  RUN  CF  FILE  STORAGE 


COWON/SUR  F/PSUR  ( 161) 

C0M4ON/FIX/0MEGA 
COMMON/MU  TUR  /CMU  (161.41) 

C0U40N/SK I NCF/CF( 1 6 1 ) 

C0M40N/GRID1/X(161 .41) ,Z(161 ,41) 

COfcMON/PAR/GAMHA , REYREF .ALFA, ALFAI .REDFRE . AMINF .ALFAI 
C0M4GN/DGRID/DT. IMAX, KMAX, ITEL,  ITEU 
C0M40N/GR I O/YACOB (161.41) 

COfcMON /DAMP /WW , WW I 

C0M40N/FL0W/Q1 (161 .41) .02(161 .41) .03(161 .41). 04(1 61 ,41) 
DIMENSION  TITLE (15) ,TYTLE(15) . ALPHA(96 ) ,CPTH( 97 , 96) . XTH(97) 
DIMENSION  XR( 161 ,49)  ,ZR(161 , 49) . 01R( 1 61 . 41 ) , Q2R( 1 61 . 41 ) 
DIMENSION  Q3R( 161 .41 ) , 04R( 161 .41 ) 

C0M4ON/LOC  I C/RSTRT  .PITCH,  PLUNGE 
LOGICAL  RSTRT. PITCH. PLUNGE 

READ  INPUT  DATA 


READ  (5.5) 

READ  (5,1)  TITLE 
READ  (5.5) 

READ  (5.2)  IMAX. KMAX, DT.WW, ALFA, ALFAI .ALFAI .REDFRE, AMINF 
READ  (5.5) 

READ  (5.3)  ITEL. ITEU. REYREF. DNMIN.TSTART. FORMAT. RSTRT, PITCH. PLUNGE 
READ  (5.5) 

READ  (5.4)  CSTP.CPLT.NSTP.PSTP 

INITIALIZE 

GAMMA  —  1.4 
PI  «  ATAN( 1 . )«4. 

SET  ALFAD  FOR  STEADY  STATE  PL0T3D  OUTPUT 
ALFAD  -  ALFA 

FORCE  DT  TO  BE  EQUAL  TO  UNITY  FOR  STEADY  FLOW  PROBLEMS 
THIS  INVOKES  THE  RELAXATION  MOOE 
IF(REDFRE. LE. 0.001)  DT  -  1.0 
REYREF  -  REYREF  •  1 . E+06 
NITN  -  CSTP  +  NSTP 
NPLOTS  -  CPLT 
CPLT  -  CPLT  +  1 
CSTP  -  CSTP  +  1 

DENSITY  NORMALISED  WITH  RESPECT  TO  ROINF 
VELOCITIES  NORMALISED  WITH  RESPECT  TO  AINF 
TOTAL  ENERGY  NORMALISED  WITH  RESPECT  TO  (ROINF.AINF.AINF) 
TOTEN-AMINF. AMINF. 0,5+1 ./( GAMMA •( GAMMA- 1 . )) 

ALFA  -  ALFA  .  ATAN(1.)  /  45. 

ALMEAN  -  ALFA 

ALFAI  -  ALFAI  •  ATAN( 1 . )  /  45. 

ALFAI  -  ALFAI  «  ATAN( 1 . )  /  45. 

ALFAS  «  ALMEAN  -  ALFAI 

WR I TE(6 )  '  PLOT  ITN  TIME  AOA  CL  CD  CM' 

UINF  -  AMINF 
VINF  -  0.0 
DO  7  1-1 . IMAX 
DO  7  K-l , KMAX 
01(1 , K )— 1 . 

02 ( I . K)-UINF 
03 ( I . K )— VI NF 
04(1 ,K)-TOTEN 
CONTINUE 

INPUT  STORED  FLOW  FIELD  DATA 

IF(RSTRT)THEN 

I F ( FORMAT .  LT . 0 . 0  ) THEN 

READ( 7)  TIME, 01 ,02,03.04 

ELSE 

READ( 7 )  IMAX, KMAX 


ooo  -  oooo 


READ(7)  AMINF.ALFAD.REYREF.TIME 
READ(7)  ({01(1, J).  I  — 1  . I  MAX) ,  J-1  . KMAX ) , 

+  ((02(1. J).  1-1 . IMAX) ,J-1 .KMAX) , 

+  ((03(1, J).  1-1. IMAX),J-1. KMAX), 

+  ((04(1. J).  1-1 . IMAX) .J-1 .KMAX) 

END  IF 

IF(NPLOTS.GT ,0)THEN 
DO  9  K1-1.NPL0TS 

READ(90)  IMAXR , KMAXR 

READ(90)  ((XR(I.J),  J-1. IMAXR). J-1. KMAXR), 

+  ((ZR(I.J),  1-1 . IMAXR). J-1 .KMAXR) 

WRITE(20)  IMAXR, KMAXR 

WRITE(20)  ((XR(I.J),  1-1 . IMAXR) .J-1 .KMAXR) . 

+  ((ZR(l.J),  1-1 . IMAXR) , J-1 .KMAXR) 

9  CONTINUE 

DO  10  K2-1.NPL0TS 

READ(91 )  IMAXR , KMAXR 

READ(91 )  AMINFR.ALFADR.REYREFR.TIMER 

READ(91 )  ((QIR(I.J).  1-1 . IMAXR) . J-1 .KMAXR) , 

+  ( (Q2R( I . J ) .  1-1 , IMAXR) .J-1 .KMAXR) , 

+  ( (03R( I . J ) ,  1-1 . IMAXR) .J-1 .KMAXR) , 

+  ( (Q4R( I , J ) ,  1-1 . IMAXR) . J-1 .KMAXR) 

WRITE(21 )  IMAXR, KMAXR 

WRIT£(21)  AMINFR.ALFADR.REYREFR.TIMER 

WRITE(21)  ((QIR(I.J),  1-1. IMAXR), J-1 . KMAXR), 
+  ( (02R( I . J) .  1-1 . IMAXR) . J-1 .KMAXR) . 

+  ( (Q3R( I . J )  .  1-1. IMAXR), J-1, KMAXR). 

+  ( (04R( I . J ) ,  1-1 . IMAXR) .J-1 .KMAXR) 

10  CONTINUE 

DO  11  K3-1.NPL0TS 

READ (92)  IMAXR. KMAXR 

READ(92)  ((XR(l.J).  1-1 .IMAXR) , J-1 .KMAXR) . 

+  ((ZR(l.J),  1-1 .IMAXR) .J-1 .KMAXR) 

WRITE(50)  IMAXR. KMAXR 

WRITE(50)  ((XR(I.J).  1-1, IMAXR). J-1 .KMAXR), 

+  ((ZR(I.J).  1-1, IMAXR). J-1 .KMAXR) 

11  CONTINUE 

DO  12  K4-1 .NPLOTS 

READ (93)  IMAXR. KMAXR 

READ(93)  ((XR(I.J).  1-1 . IMAXR) . J-1 . KMAXR) , 

+  ((ZR(I.J).  1-1 .IMAXR). J-1 .KMAXR) 

WR I T  E ( 60 )  IMAXR. KMAXR 

WR I T E ( 60 )  ((XR(I.J),  1-1 .IMAXR) .J-1 .KMAXR) , 

+  ((ZR(I.J).  1-1 .IMAXR), J-1 .KMAXR) 

12  CONTINUE 
END  IF 

END  IF 

I F( TSTART  GE . 0 ■ )  TIME  -  TSTART 
I F( . NOT . (RSTRT) )  TIME  -  0. 

•••  GENERATE  COMPUTATIONAL  GRID.  DEFINE  SURFACE  COORD  •  0  AOA, 
ROTATE  TO  INITIAL  AOA  AND  COMPUTE  METRICS 

CALL  AIRFOL( IMAX. KMAX. ITEL, ITEU) 

IF  (REYREF.GT.0)  CALL  CLUSTR(DNMIN) 

TEDGE  -  X ( ITEL  +  48,1) 

DO  15  I  -  1.97 

XTH(I)  -  X( I+ITEL-1 . 1 ) -TEDGE 
5  CONTINUE 

CALL  ROTGRID(X.Z, IMAX. KMAX. ALFAS) 

CALL  METRIC 

...  ITERATIVE  LOOP 


DO  1000  ITN-1.NSTP 
TIME  -  TIME  +  DT 


ooo  o  cn  cn  >  o<  ro  —  O  -•  OOO  OOO  OOO  OOO  oo 


...  ROTATE  GRID  TO  NEW  ANGLE  OR  SET  TO  CORRECT  ANGLE  FOR  RESTARTS 
IF  (PITCH)  THEN 

OMEGA  «  2.  •  REDFRE  •AMINF»SIN(REDFRE  •  2.  *  TIME  •  AMINF) 
1  «ALFA1 

ALOLD  -  ALMEAN  -  ALFA1  •  C0S(2.  •  REDFRE  •  AMINF  • 

1  (TIME  -  DT)) 

ALFA  -  ALMEAN  -  ALFA1  •  COS(REDFRE  •  2.  •  TIME  •  AMINF) 
ALFAD  -  ALFA  •  45 .  /  ATAN( 1 . ) 

DALFA  -  ALFA  -  ALOLD 

IF(RSTRT.AND.ITN.EQ.I)  DALFA  »  ALFA  -  2..ALFAS 
CALL  R0TGR1D(X.Z. IMAX.KMAX. DALFA) 

CALL  METRIC 
END  IF 

IF  (PLUNGE)  THEN 

OMEGA  -  2.  •  REDFRE  •  AMINF 
ALFA  -  ALMEAN 
END  IF 

...  COMPUTE  THE  SOLUTION  BY  ADI  TECHNIQUE 
CALL  SLPS( I TN, OMEGA, DALFA) 

...  APPLY  BOUNDARY  CONDITIONS 
CALL  WALLBC 

...  GENERATE  PLOT3D  FILES 

IF(CSTP.EO. (CPLT.PSTP))  THEN 
WRITE(20)  IMAX.  KMAX 

WRITE(20)  {(X(I.J).  1=1 . IMAX ) ,  J-1.KMAX), 

1  ((Z(I.J).  1=1, IMAX).  J-1.KMAX) 

WRITE(21)  IMAX.  KMAX 

WRITE(21)  AMINF.  ALFAD,  REYREF,  TIME 

WRITE(21 )  ((01(1. J).  1=1. IMAX).  J-1.KMAX), 

1  ((02(1. J).  1=1. IMAX).  J=1 . KMAX ) , 

2  ( ( 03 ( I . J ) .  1=1. IMAX),  J-1.KMAX). 

3  ((04(1. J),  1=1, IMAX),  J=1,KMAX) 

...  GENERATE  PERFORMANCE  COEFFICIENTS 

CALL  LOAD(PSUR.CL.CD.CM.ALFAS) 

AOA  =  ALFA. 180. /PI 

WRITE! 6. 6)  CPLT .CSTP. TIME. AOA, CL. CD, CM 

CALL  CPPLOT (ITEL.ITEU.AMINF.X(I.I).Z(I .1),PSUR) 

CPLT  =  CPLT  +  1 
END  IF 

CSTP  .  CSTP  +  1 
300  CONTINUE 

FORMAT  (1X.15A4) 

FORMAT  (1X.2I7.7F8.4) 

FORMAT  ( 1 X . 2 1 7 , 4F8. 5 ,3L7) 

FORMAT  (IX, 417) 

FORMAT  (IX) 

FORMAT  (IX, 13. I6.F8.3.F9.5.3F8.4) 

END 


SUBROUTINE  AMAT1 (K, IMX1 .XIX.XIZ.XIT) 

COAWON/F LOW/01 (161,41). Q2 (161 .41). 03(161 ,41 ) ,04( 1 61 . 41 ) 
C0AW0N/PERTR/DQ1 (161,41), DQ2( 1 61 , 41 ) , DQ3( 1 6 1 , 41 ) , DQ4( 1 61 , 41 ) 
CCA440N/ AM/ A  (  4 . 4 , 1  61  ) 

COAMON/PAR/GAM.IA , REYREF  ,  ALFA .  ALFA1  .REDFRE  ,  AMI NF . ALFAI 


ooo  ooo  ooo 


DIMENSION  XIT( 161 .41),XIX(161,41),XIZ(161,41) 

REAL  K0.K1 . K2 

...  AMAT1  COMPUTES  THE  COEFFICIENT  MATRIX  DE/DO  DURING  XI  SWEEP 

GUI  -  GAM4A  -  1 . 

DO  1000  I  -  2  .  IMX1 

K0  -  X I T ( I , K) 

K1  -  XIX( I ,K) 

K2  -  XIZ(I.K) 

U  -  02(1, K)  /  01(1. K) 

W  -  03(1. K)  /  01(1. K) 

EBYR  -  04 ( I . K )  /  01(1 . K) 

PHI 2  -  0.5  .  GM1  *  (U  •  U  +  W  •  W) 

THETA  -  K1  •  U  +  K2  «  W 

A(1 .1 .1)  -  K0 

A(1 .2.1)  -  K1 

A(1 ,3.1)  «  K2 

A(1 .4. 1)  -  0 

A(2 . 1 . I )  -  K1  •  PH 1 2  -  U  »  THETA 

A(2 . 2  -  K0  +  THETA  -  K 1  .  (GM1  -  1.)  •  U 

A(2 . 3  K2  *  U  -  GM1  *  K1  •  W 

A(2.«  <1  *  CM1 

A(3 . '  <2  *  PHI2  -  W  •  THETA 

A(3.:  <1  *  W  -  K2  *  GM1  «  U 

A(3.  ■  *0  +  THETA  -  K2  *  (GUI  -  1 . )  «  W 

A(3 .  -  *2  *  GM1 

A(4 , :  «  THETA  .  (2.  .  PHI2  -  GAMMA  •  EBYR) 

A(4.:  .  K1  •  (GAM4A  .  EBYR  -  PHI2)  -  GM1  .  U  •  THETA 

A(4 , 3 .  -  K2  .  (GAM4A  .  EBYR  -  PHI2)  -  GM1  .  W  .  THETA 

A(4.4..  -  K0  +  GAM4A  *  THETA 

1000  CONTINUE 
RETURN 
END 


SUBROUTINE  AMAT2( I ,KMX1 , ZETAX . ZETAZ . ZETAT) 

C0M40N/F LOW/01 (161, 41), 02 (161 . 41 ) ,03( 1 61 . 41 ) ,04( 1 61 , 41 ) 
C0M40N/PERTR/DQ1 (161 .41). 002(161 .41). 003(161 ,41).D04(161 .41) 
C0M40N/ AM/ A  ( 4 , 4 . 1 6 1  ) 

C0M40N/PAR/GAIH4A .  REYREF ,  ALFA .  ALFA1  .  REDFRE ,  AMI  NF ,  ALFAI 
DIMENSION  ZETAX (161, 41). ZETAZ (161 .41 ) .ZETAT( 161 ,41) 

REAL  K0.K1 , K2 

..  AMAT2  COMPUTES  THE  COEFFICIENT  MATRIX  DF/DQ  DURING  ETA  SWEEP 

GM1  -  GAX44A  -  1  . 

DO  1000  K  *  2  ,  KMX1 
K0  -  ZETAT(I.K) 

K 1  -  ZETAX(I.K) 

K2  -  ZETAZ(I.K) 

U  -  Q2 ( I . K )  /  01(1 ,K) 

W  -  03(1, K)  /  01(1. K) 

EBYR  -  04(1. K)  /  01(1. K) 

PH12  -  0.5  •  GM1  .  (U  •  U  +  W  »  W) 

THETA  -  K1  .  U  +  K2  *  W 
A( 1 . 1 ,K)  -  K0 

A( 1 , 2 ,K)  -  K1 

A(1.3.K)  -  K2 

A(1.4.K)  -  0 

A( 2 . 1 ,K )  -  K1  .  PH! 2  -U  «  THETA 

A(2.2.K)  -  K0  +  THETA  -  K 1  •  (GM1-1 . )  *  U 

A(2.3,K)  -  K2  »  U  -  GM1  .  K1  .  W 

A( 2 . 4 , K )  -  K1  •  GM1 

A(3.1.K)  -  K2  «  PHI2  -  W  .  THETA 

A(3.2.K)  -  K1  •  W  -  K2  *  GM1  •  U 

A(3 . 3  ,K)  -  K0  -*■  THETA  -K2  .  (GM1-1  .  )  •  W 
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A(3 . 4 ,K) 
A(4 . 1 .K) 
A(4,2.K) 
A( 4 , 3 , K) 
A( 4 , 4 , K ) 

1000  CONTINUE 
RETURN 
END 


-  K2  •  GUI 

-  THETA  .  (2.  •  PHI2  -  GAMMA  •  EBYR) 

-  K1  .  (GAMMA  •  EBYR  -  PHI2)  -  GM1  .  U  •  THETA 

-  K2  *  (GAMMA  •  EBYR  -  PHI2)  -  GM1  .  W  .  THETA 

-  K0  +  GAMMA  •  THETA 


SUBROUTINE  SLPS( I TN . OMEGA , DALFA) 

COMMON/ FL0W/Q1 ( 161 .41 ) .02(161 .41). Q3( 161 .41 ) .04(161 ,41) 
C0MM0N/P£RTR/DQ1( 161 . 41  )  .  DQ2(  161  ,41  )  .003(161  ,41  )  ,004(161  .41) 
COMMON/ AM/ A ( 4 . 4 , 161) 

COMMON/TRID/DD( 4 .4,161 . 41 ) ,MM(4 . 4 . 1 61 .41 ) . EE(4.4. 161 .41) 

1 . GG( 4,161.41) 

COMMON/PAR/GAMMA  .REYREF . ALFA . ALFA1 . REDFRE . AMI NF . ALFAI 
COMMON/DGR I D/DT . I MAX , KMAX . ITEL, ITEU 
COMMON/GR I D/YACOB (161.41) 

COMMON/DAMP/WW ,  WW I 

COMMON/MTR1X/  XIX (161 .41 ) . X I Z( 1 61 . 41 ) . 2ETAX ( 1 6 1 .41) . 
+ZETAZ(161.41) 

1  , X  I T ( 1 61 ,41).ZETAT(161 .41) 

REAL  MmI 

DIMENSION  DELTAT(161 .41) 

SUBROUTINE  SLPS  DOES  THE  BULK  OF  THE  WORK  FOR  THE  AC!  ALGORITHM. 
IT  CALLS  FLUX  AND  COMPUTES  RIGHT  HAND  SIDE  DURING  THE  TWO  SWEEPS, 
ASSEMBLES  THE  COEFFICIENT  MARICES.  ADDS  IMPLICIT  AND  EXPLICIT 
DISSIPATION  AND  CALLS  THE  TRIDIAGONAL  SOLVER  TO  OBTAIN  THE  FINAL 
SOLUTION. 

!  ! SET  VALUE  OF  IMPLICIT  DAMPING  COEFFICIENT 
WW1  «  20.0  •  WW 

IM1  -  IMAX  -  1 

IM2  -  IMAX  -  2 

KM1  «  KMAX  -  1 

KM2  -  KMAX  -  2 

IF(ITN.EO.I)  THEN 

!! LOCAL  TIME  STEP  OPTION  FOR  STEADY  STATE  CALCULATIONS 
1F(REDFRE.LT. 0.001)  THEN 
DO  777  K  -  2  ,  KMAX  -  1 

DO  777  I  »  2  .  IMAX  -  1 


DO  777  K  -  2 
DO  777  I  »  2 
DELTAT(I.K)  - 
CONTINUE 
ELSE 

DO  778  K  «  2 
DO  778  I  -  2 
DELTAT ( I , K)  » 
END  IF 
END  IF 


0.5  /  (1.  +  SORT(ABS(YACOB(I ,K)))) 


,  KMAX  -  1 
IMAX  -  1 
1 .0 


•••  THE  DISSIPATION  TERMS  ARE  CONSTRUCTED  AND  STORED  IN  THE  ARRAYS  DQ1 . 
D02.DQ3  AND  004. 

CALL  DISSIP 

...  THE  RIGHT  HAND  SIDE  AT  KNOWN  TIME  LEVEL  IS  NOW  COMPUTED  AND  ADDED 
CALL  RESI (OMEGA) 

...  IF  VISCOUS  FLOW  IS  COMPUTED  THE  VISCOUS  TERMS  ARE  ADDED  TO  DQ1  ETC.  HERE. 

I Ff REYREF . GT . 0 . )  CALL  STRESS( ITN.DALFA) 

...  I -SWEEP. 

DTH  -  DT  .  0.5 
DTW  -  DT  .  WWI 


$ 


f-t  f,*  j.t  , 


DO  3  K  -  2  .  KM1 

CALL  AMAT 1 (K, I MAX- 1 ,XIX,XIZ,XIT) 
DO  4  11  -1.4 

DO  4  12  -1.4 

00  5  I  -  2  .  I MAX  -  1 

£E{!1 .12.1-1 ,K)  -  A( 11.12. 1  +  1 ) 

DD(I1 .12.1-1 .K)  -  -  A( 1 1 . 12. 1—1 ) 
5  CONTINUE 
4  CONTINUE 

•  IMPLICIT  DAMPING  ADDED  HERE. 


DTK  .  DELTAT(I.K) 
DTH  .  DELTAT ( 1 tK) 


DO  6  II 
DO  6  I 
DTI 

DO (11  .11.1-1 .K) 
£E(I1.:i.I-1.K) 
1*4(11  .  II.  1-1  .K) 
6  CONTINUE 
DO  990  I 
GG( 1 , 1-1. K) 
GG(2. 1-1 .K) 
GG(3. I-1.K) 
GG(4. I-1.K) 

990  CONTINUE 
3  CONTINUE 


-1.4 

-  2  ,  1MAX  -  1 

-  DTX  /  YACOB(I.K)  •  OELTAT(I.K) 

-  DD( 11.11. 1—1 .K)  -  DTI  .  YACOB(I-I.K) 

-  ££(11,11, 1-1, K)  -  DTI  .  YACOB( 1+1 ,  K) 

-  1.  +  2.  •  OTW  .  DELTAT ( I ,K) 

-  2  I MAX  -  1 

-  DOI(I.K)  •  DELTAT ( 1 , K) 

-  D02( I ,K)  •  DELTAT ( I ,K) 

-  DQ3 ( 1 . K )  •  DELTAT ( I .K) 

-  D04(I.K)  .  DELTAT ( I , K) 


PERFORM  BLOCK  TRIDIAGONAL  MATRIX  INVERSION  FOR  THE  ENTIRE  PLANE 

CALL  MATRX 1 ( I MAX . KMAX ) 

DO  991  K  -  2  ,  KMAX  -  1 

DO  991  I  -  2  .  IM1 

DQI(l.K)  -  GG(I.I-I.K) 

D02 ( I . K )  -  GG(2, 1-1 .K) 

DQ3(I.K)  -  GG(3 . 1-1 . K ) 

D04 ( I . K )  -  GG(4, 1-1 ,K) 

991  CONTINUE 

K-SWEEP  BEGINS  HERE. 

DO  13  I  -  2  .  IM1 

CALL  AMAT2( I , KMAX— 1 .ZETAX . 2ETAZ . ZETAT) 

DO  15  II  -1,4 

DO  15  12  -1.4 

DO  15  K  -  2  .  KMAX  -  1 

EE( 1 1 .12.1 .K-1)  -  A ( 1 1  ,12, K+1 ) *DTH  •  DELTAT(I.K) 

D0(  11.12,1 ,K-1 )  -  — A( 1 1 .12. K-1). DTH  .  DELTAT (I ,K) 

15  CONTINUE 

SECONO  ORDER  DAMPING  ADDED  HERE. 


DO  16  II  -1.4 

DO  16  K  -  2  .  KMAX  -  1 

DTI  -  DTR  /  YACOB(I.K)  •  DELTAT(I.K) 

DD( 1 1.1 1.1, K-1)  -  00(11, II. I, K-1)  -  DTI  »  YACOB(I.K-l) 
EE(I1 .11. I. K-1)  -  EE(I1 . II . I .K-1 )  -  DTI  •  YACOB( I ,K+1 ) 
16  1*4(11. II. I. K-1)  -  1.  +  2.  •  DTKI  •  DELTAT(I.K) 

DO  17  K  -  2  ,  KMAX  -  1 

GG ( 1 . I .K-1)  -  DQ1 ( I . K ) 

GG( 2 . I , K— 1 )  -  DQ2 ( I . K j 

GG(3 . I ,K— 1 )  -  DQ3(I.K) 

GG(4 , 1  ,K— 1 )  -  D04(I.K) 


DO  17  K 
GG ( 1  .  I  .K-1) 
GG ( 2 . I .K-1) 
GG(3.  I. K-1) 
GG ( 4 , 1 .K-1) 
17  CONTINUE 
13  CONTINUE 


PERFORM  BLOCK  TRIDIAGONAL  MATRIX  INVERSION  FOR  THE  ENTIRE  PLANE 


I 

£ 

i 


.‘V.vv.  .'.  ^ y.  .\y  v-,-. 


W.v, 


oooo  ooo  noooo  o  oo 


CALL  MATRX2( IMAX.KMAX) 


00  18  I 
DO  18  K 
001(1 ,K) 
002(1 ,K) 
003(1 .K) 
D04(I.K) 
18  CONTINUE 


-  2  .  I MAX  -  1 

-  2  .  KM1 

-  GG(I.I.K-I) 

-  GG(2, I . K— 1 ) 

-  GG( 3 , I . K— 1 ) 

-  GC(4,I,K-1) 


YACOB(I.K) 
YAC0B( I .K) 
YAC0B( I .K) 
YAC0B( I .K) 


•••  UPDATE  FLOW  VARIABLES  AT  INTERIOR  POINTS. 

967  CONTINUE 

RMAX  -  0 . 

RUMAX  -  0 . 

RVMAX  -  0  . 

EMAX  *=  0  . 

00  995  K  -  2  .  KM1 

DO  19  I  -  2  ,  IM1 

01(1. K)  -  01(1. K)  +  001(1. K)  •  YACOB(I.K) 

02(1. K)  -  02(1, K)  +  D02(I.K)  •  YACOB(I.K) 

03(1. K)  -  03(1. K)  +  DQ3 ( I , K )  •  YACOB(I.K) 

04(1. K)  -  04(1, K)  +  004(1. K)  .  YACOB(I.K) 

19  CONTINUE 

I  !M  (DETERMINE  WHERE  IN  FLOW  FIELD  DENSITY  IS  CHANGING  MOST  RAPIDLY 
DO  995  I  -  2  .  IMAX  -  1 

IF  (RMAX . LT . ABS(DQ1 ( I , K ) «YACOB( I . K) ) )  THEN 
IR  -  I 

KR  -  K 

END  IF 

RMAX  -  AMAX1(RMAX,ABS(D01(I.K)  •  YACOB(I.K))) 

RUMAX  -  AMAX1 ( RUMAX , ABS (DQ2 ( I ,K)  •  YAC08( I . K ) ) ) 

RVMAX  -  AMAX1  (RVMAX. ABS (003(1 .K)  .  YACOB(I.K))) 

EMAX  «  AMAX1 ( EMAX. ABS (DQ4(! ,K)  •  YACOB(I.K))) 

995  CONTINUE 

I F( ( ITN-1 )/ 100* 100. EO. ( I TN— 1 ) )  WRITE  (6,3002) 

I F( ITN  EO.  0)  WRITE  (6,3002) 

!  !  !  ! ! SELECT  INTERVAL  AT  WHICH  OUTPUT  OF  RESIDUALS  IS  DESIRED 

I  F(  ( ITN-1 )/10*10.EQ.( ITN-1 ) )  WRITE  (6.3001)  RMAX . RUMAX . RVMAX , 

1  EMAX  .  IR  .KR 
RETURN 

5002  FORMAT (//. 4X , ' DRMAX ‘ ,  1 1 X . 'OUMAX' . 1 1X. 'DVMAX' , 1 1X, 'OEMAX'.10X, 

1 1 IR  ’  ,3X , ‘ KR ’ ) 

5001  FORMAT (4(E14.8.2X),2I5) 

END 


ij'  i.1  _i 


XWTOJWWoVL'W^VOTkV.WJ 


oo  ieee  I  -  2  .  imax  -  2 
do  5  n-1,4 

DO  2  K  -  2  .  KMAX  -  1 
C ( 1 1  . 1  ,K)  -  GG(II.I.K)  -  00(11.1. I. K) 


2  CONTINUE 

DO  5  12  -  1  .  4 

DO  5  K  -  2  .  KMAX  -  1 

A ( 1 1  ,  I  2  ,K)—  144(11 ,12.1. K)  -  DD(  1 1  . 

1  -  00(11  , 

2  -  DD(I1  . 

3  —  DD( 1 1 . 

C ( 1 1 . 12+1 ,K)«  EE(I1 .I2.I.K) 

5  CONTINUE 

DO  3  K  -  2  .  KMAX  -  1 
L1 1  -  A(1 . 1 .K) 

LI  I  —  1 .  /  L1 1 
U12  -  A(1 ,2.K)  •  LI  I 
U13  -  A ( 1  ,  3  ,  K  )  •  LI  I 
U1*  «  A(1 ,4,K)  •  LI  I 
L21  -  A(2. 1  . K ) 

L31  -  A(3, 1 .K) 

L41  -  A(4 , 1 ,K) 

L22  -  A(2 , 2 , K)  -  L21  •  U12 
L2I  -  1 .  /  L22 

U23  -  (A(2.3.K)  -  L21  •  U13)  •  L2I 

U24  -  ( A( 2 . 4 , K )  -  L21  •  U14)  .  L2I 

L32  -  A(3.2.K)  -  L31  *  U12 

L42  -  A(4. 2 , K)  -  L41  .  U12 

L33  -  A(3.3.K)  -  L31  •  U13  -  L32  • 

L31  -  1 .  /  L33 

U34  -  ( A(3 . 4 , K)  -  L31  •  U14  -  L32 


-  DD( 1 1 . 1 . I ,K)  .  GG(  1  , 1-1 , 

-  00(11.2.1  , K )  .  GG(2, 1-1 , 

-  00 (  II .3.  I  .K)  •  GG(3. 1-1 , 

-  00(11.4.1 , K )  •  GG(4. 1-1 , 


-  00(11.1. I, K)  •  HH( 1.12, 

-  DD( 1 1 .2. I ,K)  •  HH(2 . 12. 

-  DD( 1 1 ,3. I ,K)  •  HH ( 3 , 12. 

-  DD( 1 1 .4. I  ,K)  •  HH(4, 12. 


1-1  .K) 
1-1 .K) 
1-1 .K) 
1-1  .K) 


U14)  .  L2I 
U12 
U12 

U13  -  L32  •  U23 


U34  -  ( A(3 , 4 . K)  -  L31  •  U14  -  L32  •  U24)  •  L3I 

L43  -  A(4.3.K)  -  L41  •  U13  -  L42  •  U23 

L44  -  A(4 , 4 , K)  -  L41  .  U14  -  L42  •  U24  -  L43  •  U34 

L41  —  1 .  /  L44 

C(1.1.K)  -  C ( 1  .1 . K )  .  Lit 

C(2.1.K)  -  (C(2.1.K)  -  L21  .  C( 1 , 1 ,K>)  «  L2I 

C(3.1.K)  -  (C(3. 1 ,K)  -  L31  •  C(1.1.K) 

I  -  L32  .  C(2 . 1 , K ) )  .  L3I 

C ( 4  1.K)  -  (C ( 4 . 1 . K )  -  L41  .  C(I.I.K)  -  L42  •  C(2.1, 
I  -  L43  .  C(3 , 1 . K) )  . 

C  - . K )  -  C(1 . 2 . K )  •  L1I 

C  K)  -  (C( 2 .2 . K)  -  L21  .  C( 1 . 2 ,K) )  •  L2I 

C  .K)  -  (C(3 .2 , K)  -  L31  •  C( 1 . 2 , K) 

t  -  L32  •  C(2 . 2 , K) )  .  L3I 

C (4 , 2 ,K)  -  (C(4.2.K)  -  L41  «  C( 1 . 2 . K)  -  L42  •  C(2,2. 

-  L43  •  C (3 , 2 , K ) )  . 

C(1.3.K)  -  C( 1 . 3 , K)  .  LI  I 

C ( 2 .3 . K)  -  (C( 2 . 3 . K)  -  L21  •  C(1.3,K))  .  L2I 

C(3.3.K)  -  (C(3 , 3 , K)  -  L31  «  C(1.3.K) 

-  L32  .  C(2 . 3 , K) )  .  L3I 

C ( 4 , 3 , K)  -  (C(4,3,K)  -  L41  «  C( 1 . 3 , K)  -  L42  •  C(2,3. 

-  L43  .  C(3 , 3 , K ) )  . 

C ( 1 ,4,K)  -  C(1 .4.K)  .  LI  I 

C(2.4.K)  -  (C ( 2 . 4 . K )  -  L21  •  C ( 1 . 4 , K)  )  •  L2I 

C(3,4,K)  -  (C(3.4 . K)  -  L31  •  C(1.4,K) 

-  L32  •  C( 2 , 4 , K) )  •  L3I 

C(4,4,K)  -  (C( 4 , 4 . K )  -  L41  .  C(1,4.K)  -  L42  •  C(2.4. 

-  L43  •  C(3 , 4.K) )  . 

C(1 . 5 . K)  -  C(1 .5.K)  .  LI  I 

C(2 . 5 , K )  -  (C(2.5.K)  -  L21  •  C(1.5.K)>  .  L2! 

C( 3.5. K)  -  (C( 3 . 5 , K )  -  L31  .  C(1,5,K) 

-  L32  •  C(2 . 5 , K)  )  .  L3I 

C(4.5.K)  -  (C(4 . 5 . K)  -  L41  .  C( 1 . 5 , K )  -  L42  •  C(2,5. 

-  L43  .  C( 3 . 5 , K ) )  . 


C(4.4.K) 

1 

C(1 .5.K) 
C(2 ,5 . K) 
C  ( 3 . 5 .  K ) 


(C( 4 , 4 . K )  -  L41 

C(1.5,K)  .  LI  I 
(C( 2.5.K)  -  L21 
(C( 3 . 5 , K )  -  L31 


C(4.5.K)  -  (C(4 . 5 . K)  - 


ft 


6 


oooo  ooo  ooo 


C(3 . 1 , K ) 

C(3 , 1 ,K) 

-  U34 

• 

C(4 , 1 ,K) 

C(2 , 1 , K) 

0(2.1 , K) 

-  U24 

• 

C( 4 . 1 , K ) 

1 

-  U23 

• 

C(3. 1 .K) 

C(1 .1 ,K) 

C(1 ,1 ,K) 

-  U14 

• 

C( 4 . 1 ,K) 

1 

-  U13 

• 

C(3, 1 , K ) 

C(3 . 2.K) 

C(3.2.K) 

-  U3* 

» 

C(4.2.K) 

C(2.2,K) 

C(2.2.K) 

-  U24 

• 

C( 4 . 2 . K ) 

1 

-  U23 

• 

C( 3 , 2 , K) 

C(1 ,2,K) 

C(1  ,2.K) 

-  U14 

• 

C( 4 , 2 ,K) 

1 

-  U13 

• 

C( 3 , 2 ,K) 

C(3,3.K) 

C(3,3,K) 

-  U34 

• 

C( 4, 3 . K) 

C(2 .3 , K) 

C(2 .3 . K) 

-  U24 

• 

C( 4 , 3 . K) 

1 

-  U23 

• 

C( 3 , 3 , K) 

C(1 .3.K) 

C(1 . 3 , K) 

-  U14 

• 

C( 4 , 3 . K) 

1 

-  U13 

• 

C( 3 , 3 , K) 

C(3 . 4 , K) 

C(3, 4, K) 

-  U34 

• 

C  (  4 . 4 ,  K  ) 

C(2 , 4 , K) 

C(2 ,4 , K) 

-  U24 

• 

C (  4 . 4 ,K ) 

1 

-  U23 

• 

C( 3 , 4 , K) 

C(1 , 4 , K) 

C(1 . 4 . K) 

-  U14 

• 

C(4 . 4 , K) 

1 

-  U13 

• 

C( 3 , 4 ,K) 

C( 3 , 5 . K) 

C(3 . 5 , K) 

-  U34 

• 

C( 4 , 5 ,K) 

C(2.5.K) 

« 

C(2 . 5 , K) 

-  U24 

• 

C( 4 , 5 , K) 

1 

-  U23 

• 

C( 3 , 5 , K ) 

C(1 ,5,K) 

• 

C(1 .5.K) 

-  U14 

• 

C( 4 , 5 , K) 

1 

3  CONTINUE 

-  U1  3 

• 

C( 3 , 5 , K) 

00  6  11 

«  1  . 

4 

DO  9  K 

-  2  , 

KMAX  - 

1 

9  GG(II.I.K) 

-  C(  1 1 

.I.K) 

DO  6  12 

-  1  . 

4 

DO  6  K 

-  2  . 

KMAX  - 

1 

HH (11. 12. 

I.K)  -  C( I 1 

,12+I.K) 

6  CONTINUE 
1000  CONTINUE 


U12  •  C(2 . 1 ,K) 


U12  •  C(2 , 2 , K) 


U12  .  C(2 ,3,K) 


U12  •  C(2 , 4 , K) 


U12  •  C(2 , 5 , K) 


BACKWARD  SUBSTITUTION 


DO  7  I 
DO  7  II 
DO  7  K 
GG(II.I.K) 

1 

2 

3 

7  CONTINUE 
RETURN 
END 


-  I MAX  -  3,  1  .  -  1 
-1,4 

m  2  ,  KMAX  -  1 

-  GG(II.I.K)  -  HH(II.I.I.K) 

-  MH(I1.2.I,K) 

-  HH (11,3. I, K) 

-  HH (11,4, I ,K) 


GG(1 . 1+1 ,K) 
GG(2 . I+i , K ) 
GG(3, 1+1 .K) 
GG ( 4 , I + 1 , K ) 


SUBROUTINE  MATRX2(1 MAX. KMAX) 

COM40N/TRID/DD(4.4,161 . 41 ) ,k*l( 4 , 4 , 1 61 , 41 ) , EE(4 , 4 , 1  61 .41), 

1GC(4. 161 ,41 ) 

C0UA0N/SCRAT/A(4.4,161) .HH(4,4. 1 61 . 41 ) ,C(4 ,5. 161 ) 

REAL.  M4 

REAL  L1 1 .121 ,L31 ,L41 , L22, 132 . L42 . L33 , L43. L44 
2.L1I.L2I.L3I.L41 

THIS  SUBROUTINE  PERFORMS  THE  BLOCK  TRIDIAGONAL  MATRIX  IVERSION  FOR 
AN  ENTIRE  CONSTANT  PLANE  DURING  THE  ZETA-  SWEEP 

DO  1  II  -1,4 

DO  1  I  -  2  .  IMAX  -  1 

A I  -  1  /  **l(1. 1,1,1) 

GG(II.I.l)  -  GG(II.I.I)  •  AI 

HH(I1. 1.1.1)  -  EEC  1 1 . 1.1.1)  *  AI 


.V V* ,N  , 
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HH(I1 .2. 1 .1)  -  EE( 1 1 .2.1.1)  •  A I 
HH( II .3. 1 , 1 )  -  EE(  II  .3,  I  .  1  )  •  A I 
HH (11,4,1.1)  -  EE( 11,4.1,1)  «  A I 

1  CONTINUE 

DO  1000  K  -  2  ,  KMAX  -  2 
DO  5  11-1,4 

DO  2  1  -  2  .  I MAX  -  1 

C(I1.1.I)  -  GG(ll.I.K)  -  00(11,1.1 

1  -00(11,2.1 

2  -00(11,3.1 

3  -  00(11,4,1 

2  CONTINUE 

DO  5  12  -  1  .  4 

DO  5  I  -  2  .  IMAX  -  1 

A(  1 1 ,12.1)“  M4( 1 1 . 1 2 . I . K )  -  00(11  . 

1  -00(11. 

2  -  00(11  , 

3  -00(11. 

C( 1 1 .12+1 .1)-  EE( 1 1 . 12 . 1 .K) 

5  CONTINUE 

DO  3  I  -  2  .  IMAX  -  1 
L 1 1  -  A(1 . 1 . I) 


.K)  .  GG( 1 , I.K-1) 
,K)  •  GG(2. I , K— 1 ) 
.K)  .  GG(3, 1 , K-1 ) 
,K)  •  GG(4, I , K— 1 ) 


1.1. K)  •  HH( 1,12.1, K— 1 ) 

2.1. K)  •  HH(2. 12.1. K-1) 

3.1. K)  •  HH(3.I2.I.K-1) 

4. 1 . K)  •  HH(4, 12. I.K-1) 


LI  I  -  1  . 
U12  -  A( 
U13  -  A( 
U14  -  A( 
L21  -  A( 
L31  -  A( 
L41  -  A ( 
L22  -  A( 
L2 1  -  1  . 
U23  -  (A 
U24  -  (A 
L32  -  A( 
L42  -  A( 
L33  «  A( 
L3I  -  1  . 
U34  -  (A 
L43  -  A( 
L44  -  A( 
L4I  -  1  . 

C(1 .1  . 1) 
C(  2 , 1 . I) 
C(3 . 1  .  I  ) 


1  .  /  L 1 1 
A( 1 . 2  .  I  )  - 
A( 1 . 3 , I )  - 
A( 1 . 4 . I )  - 
A(2 . 1  .  I  ) 
A(3. 1 .1) 
A( 4 . 1 . I ) 
A(2 , 2 . 1)  - 
1  .  /  L22 
( A( 2 . 3 . I ) 


U1 3)  •  L2I 
U14)  .  L2 I 


( A( 2 , 4 . I )  -  L21  •  U14)  .  L2 I 
A(3 . 2 . I )  -  131  •  U12 
A( 4 , 2 . I )  -  L41  •  U12 
A(3.3.I)  -  L31  •  U13  -  L32  • 
1.  /  L33 

( A( 3 . 4 . I )  -  L31  •  U14  -  L32 
A(4.3, I)  -  L41  •  U13  -  L42  « 

A( 4 , 4 , I )  -  L41  •  U14  -  L42  • 

1  .  /  L44 

I)  -  C(1 . 1 . I )  •  LI  I 
I)  -  (C(2 . 1.1)  -  L21  •  C(1.1 

I)  -  (C(3 . 1,1)  -  L31  •  C(1,1 

-  L32  •  C(2 , 1 


•  U24)  .  L3I 
U23 

U24  -  L43  •  U34 


C( 4 , 1 . I)  -  (C(4 , 1 .1)  -  L41  .  C(1 .1 


C(1 .2.1)  -  C ( 1 . 2 . I )  •  LI  I 
C(2 . 2 .  I )  -  (C(2 . 2 . I )  -  L21 
C(3 . 2 , I )  -  (C(3,2, I )  -  L31 
1  -  L32 

C(4 , 2 . I )  -  (C(4 ,2.1)  -  L41 


C(4.2. I)  - 

1 

C(  1  . 3 . I )  - 
C( 2 . 3 , 1 )  = 
C(3 . 3 .  I )  - 

1 

C(4 . 3 , I )  - 


C(1 .3. 1 )  •  LI  1 
(C(2 . 3 , I )  -  L21 
(C(3 . 3.1)  -  L31 
-  L32 
(C(4 ,3 , I )  -  L41 


C(  1  .2 
C(1 .2 
C  ( 2 . 2 
C{  1  .2 


C(  1 .3 
C(  1  .3 
C(  2 . 3 
C(  1  .3 


C ( 1 , 4 , 1 )  -  C( 1 .4,1)  *  LI  I 

C( 2 . 4 , I )  -  (C( 2 . 4 , I )  -  L21  .  C(1 .4 

C(3 . 4  ,  I  )  -  (C(3 .4,1)  -  L31  «  C(1,4 


C(4.4, I)  -  (C(4 ,4,1)  -  L41 

1 

C(1.5,I)  -  C( 1 . 5 , I )  •  LI  I 


-  L31  «  C(1 .4 

-  L32  •  C ( 2 . 4 

-  L41  .  C(1  .4 


.!))  •  L2I 
.0 

, I))  •  L3I 

. I)  -  L42  •  C(2 .1.1) 
L43  •  C(3. 1  .  I))  •  L4I 

.0)  •  L2I 
.1) 

, I))  •  L31 

. 1)  -  L42  •  C(2.2,  I) 
L43  •  C(3 , 2 , 1 ))  •  L41 

.0)  •  L21 
.!) 

. I ) )  •  L3I 

,1)  -  L42  •  C(2,3. I ) 
L43  •  C(3 , 3 . I ))  •  L4I 

.!))  .  L2I 
•  I) 

, !))  »  L3I 

. I)  -  L42  •  C(2 , 4 , I ) 
L43  •  C( 3 , 4 , I ))  •  L4I 


C(2,5.  I) 

0(3 .5. 1 ) 


(C(2 , 5 . I ) 
(0(3,5, I ) 


C(4 . 5 , I )  -  (0(4,5, I )  - 


0(3, 1 .1) 
C(2 , 1 , I ) 


c(3.i. n 

C(2 . 1 .1) 


0(1.1.  I)  -  0(1. 1.1)  - 


0(3.2. I) 

0(2.2. I) 


0(3.2. I) 

0(2.2. I) 


0(1 .2. I)  -  0(1 .2. I)  - 


0(3.3,  I) 
0(2.3. I) 


0(3.3,  I)  - 
0(2.3. 1)  - 


0(1 .3. I)  -  0(1 .3. I)  - 


0(3.4.  I) 
0(2.4, I ) 


0(3,4. I)  - 
0(2.4. I)  - 


0(1, 4. i)  -  C( 1,4.1) 


0(3.5.  I) 
0(2,5.  I) 


0(3.5. I)  - 
0(2.5. I)  - 


0(1  .5. I)  -  0(1 ,5. I)  - 

1 

3  CONTINUE 


•  c(i,: 

•  c(i.: 

•  c(2.: 

•  c(i.: 

0(4.1, 
0(4.1. 
0(3. 1 . 
0(4.1, 
0(3.1. 
0(4,2, 
0(4.2. 
0(3.2, 
0(4,2. 
0(3.2. 
0(4,3. 
0(4,3. 
0(3,3, 
0(4,3, 
0(3.3, 
0(4.4. 
0(4.4. 
0(3.4. 
0(4,4. 
0(3.4. 
0(4,5, 
0(4,5. 
0(3.5. 
0(4.5. 
0(3.5. 


5.1))  •  12 1 
5,0 

5.1) )  •  L3I 

5.1)  -  L42  •  0(2. 5.1) 
L43  .  0(3.5. I))  •  L4I 


U12  •  0(2.1, I) 


-  U12  •  0(2.2. I) 


-  U12  •  0(2.3. I) 


-  U12  •  0(2.4, I) 


-  U12  •  0(2.5. I) 


DO  6  II 
DO  9  I 

9  GG ( 1 1 , I .K) 

DO  6  12 
DO  6  I 

HH(I1 , 12. I .K) 
6  CONTINUE 
1000  CONTINUE 


-1.4 

-  2  .  I MAX  -  1 

-  0(11.1.1) 

-1.4 

-  2  .  IMAX  -  1 

-  0(11 . 12+1 .1) 


BACKWARD  SUBSTITUTION 


DO  7  K 
DO  7  II 
DO  7  I 
GG ( 1 1 , I .K) 


7  CONTINUE 
RETURN 


KMAX  -  3.  1  .  -  1 

1  .  4 

2  IMAX  —  1 

GG(  1 1  .  I  ,K)  —  K*4(  1 1  .  1  .  I  .K)  »  GG(  1  .  I  . K+1 ) 

-  FH( 1 1 , 2 , I , K)  »  GG(2  .  I  .K+1 ) 

-  HH( II .3, I ,K)  •  GG ( 3 , I , K+1 ) 

-  HH( 1 1 .4, I ,K)  •  GG(4 , I .K+1) 


SUBROUTINE  METRIC 
COfcMON/F I X/OMEGA 

O*fct0N/DGRID/DT. IMAX, KMAX, ITEL.  ITEU 
COMMON/GR I D 1 /X (161 ,41 ) . Z ( 1 61 .41 ) 

C0M4ON/GR I D/ YACOB (161,41) 

COMMON/MTR ' X/X I  X ( 1 6 1 ,41 ) ,XIZ(161  . 41 ) . ZETAX ( 1 6 1  ,41), 

+ZETAZ( 161 .41 ) , 

1XIT( 161.41 ) ,ZETAT(161 ,41 ) 

SUBROUTINE  METRIC  COMPUTES  THE  METRICS  IN  BOTH  DIRECTIONS  AND 
THE  UNSTEADY  COEFFICIENTS  ETAT .  ETC. 
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c 

DO  1000  K  -  1  .  KMAX 
DO  1000  I  -  1  .  IMAX 
XTAU  -  OMEGA  .  Z(I ,K) 

YTAU  -  OMEGA  .  (-X( I ,K)) 

C.»»  PRESENT  SET  UP  IS  FOR  FLOW  PAST  AN  AIRFOIL. 

C 

CUM  ! CENTRAL  DIFFERENCES  AT  INTERIOR  POINTS.  TWO-POINT  ONE-SIDED 
C! ! !! (DIFFERENCES  DOWNSTREAM,  THREE-POINT  AT  OTHER  OUTER  BOUNDARIES 
I F( I  .  EQ. 1 .OR. I .£0. IMAX)  GO  TO  10 
XXI  -  .5  •  (X(I  +  1 ,  K  )  -X  ( I  —  1 , K ) ) 

ZXI  -  .5  •  (2(1+1 .K) -Z( 1-1 , K ) ) 

GO  TO  15 

10  I F ( I  .EQ. IMAX)  GO  TO  16 

XXI  -  1.0  .  (X( 2 . K )  -  X(1  , K ) ) 

ZXI  -1.0.  (Z(2,K)  -  Z(1 ,K)) 

GO  TO  15 

16  XXI  «  1.0  .  (X(IMAX.K)  -  X ( IMAX— 1 . K) ) 

ZXI  -  1.0  •  (Z(IMAX.K)  -  Z ( IMAX— 1 ,  K  )  ) 

15  CONTINUE 

I  F  ( K .  EQ .  1  .  OR .  K .  EQ . KMAX )  GO  TO  1 7 
XZET  -  .5  «(X( I .K+1 )-X( I ,K-1 ) ) 

ZZET  -  .5  - (Z( I . K+1 )-2( I . K-1 ) ) 

GO  TO  20 

17  IF(K.EO.KMAX)  GO  TO  18 

XZET  -  2.  »  X( I . 2 ) — 1 . 5  •  X(I.I)  -  .5  .  X( I .3) 

ZZET  -  2.  •  Z( I ,2)  -  1.5  •  Z(I,1)  -  .5  .  Z(I,3) 

GO  TO  20 

18  XZET  -  1.5  •  X( I . KMAX)— 2  «  X( I , KMAX-1 )+ . 5»X( I , KMAX-2) 

ZZET  -  1.5  •  Z(I ,KMAX)-2.«  Z( I .KMAX-1 )+. 5«Z( I .KMAX-2) 

20  CONTINUE 

CM!!  !  COMPUTE  JACOBIAN 

YACOBI  -  XXI  •  ZZET  -  XZET  •  ZXI 
YACOB( I .K)  -  1 .  /  YACOBI 
XIX(I.K)  -  ZZET  .  YACOB(I.K) 

XIZ(I.K)  -  -XZET  .  YACOB(l.K) 

XTAU  -  OMEGA  .  Z(I ,K) 

YTAU  -  -  OMEGA  .  X(I ,K) 

XJT(I.K)  -  -  XIX(I.K)  .  XTAU  -  XIZ(I.K)  •  YTAU 
ZETAX(I.K)  -  -ZXI  •  YACOB(I.K) 

ZETAZ(I.K)  -  XXI  .  YACOB(I.K) 

ZETAT(I.K)  -  -  ZETAX(I.K)  .  XTAU  -  ZETAZ(I.K)  .  YTAU 
1000  CONTINUE 
RETURN 
END 

C 


- - -  - - -  ~ - ’  ~  ^ - - - - ~ - -  ”  ~  - - ~  ~ 

c 

SUBROUTINE  DISSIP 

C0M4ON/FLOW/Q1 (161 ,41 ) , Q2( 1 6 1 . 41 ) . Q3( 1 6 1 . 41 ) . Q4 ( 1 61 . 41 ) 
C0M4ON/PERTR/DQ1 (161 . 41 ) ,DQ2 ( 1 61 , 41 ) , DQ3( 1 61 , 41 ) ,DQ4( 1 61 , 41 ) 
CONMON/DGR I D/DT . IMAX. KMAX. I  TEL. I TEU 
COMMON /GR I D/ Y ACOB (161.41) 

COMMON/DAMP/WW ,  WW I 

DIMENSION  P( 161 ) .EPS(T6T).DIS1 (161 , 4) . D I S2( 1 6 1 , 4) 

C 

C  this  SUBROUTINE  ADDS  THE  FOURTH  ORDER  DISSIPATION  TERMS  TO  THE 
C  RIGHT  HAND  SIDE 

C 

I  Ml  -  IMAX  -  1 

KM1  -  KMAX  -  1 

I  M2  «  IMAX  -  2 

KM2  -  KMAX  -  2 

DO  10  K— 2  .  KM1 

:  COMPUTE  SWTICHING  FUNCTION  BASED  ON  SECOND  DERIVATIVE  OF  PRESSURE 

DO  1  I  -  1  .  IMAX 


1  P(I)  .  .4  .  (Q4(I.K)-(Q2(I . K) »»2+Q3( I .K).*2)/(2..Q1( I . K ) ) ) 

DO  2  I  -1  .  I MAX 

IP2  -  I  +  2 

IF(I.EQ.IMI)  IP2  -  IMAX 
IM2  -  1  -  2 
I F  ( I .£0.2)  I M2  -  1 
IP1  -  I  +  1 

I  F( I  .EQ. IMAX)  IP1  -  IMAX 

IM  -  I  -  1 

I F( I  E0.1)  IM  -  1 

I F ( I . EO. 1 )  IM2  -  1 

IF( I . EO. IMAX)  IP2  -  IMAX 

EPS ( I )  -  ABS(P(IP1)-2..P(I)+P(IM))/ABS(P(IP1)+2..P(I)+P(IM)) 

C  VOL  -  2.  /(YACOB( I ,K)+YACOB( I P 1 ,K) ) 

VOL  -  1 . 

DISI(I.I)  -  (Q1 ( IP1 . K)— 01 ( I . K ) ) *VOL 

DIS1 ( I . 2)  -  ( 02 ( I P 1 ,K)-Q2( I ,K) ) *VOL 

DIS1 ( I .3)  -  ( 03 ( I P 1 , K)— Q3( I ,K) ) «V0L 
HP1  -  04  ( I  PI ,K)+P(IP1) 

HP  -  04 ( I , K )+P{ I ) 

HM1  -  04( IM.K)  +  P( IM) 

HP2  -  Q4( IP2 . K)  +  P(IP2) 

HPM  =  04 ( IM , K  )+P ( IM) 

DIS1 (I .4)  -  ( HP  1 —HP ) • VOL 

DIS2(I,1)  —  (01(IP2.K)-3.»(Q1(IP1 ,K)-01 (I.K))— Q1(IM,K) )«V0L 
DIS2( I . 2)  —  (Q2(IP2,K)-3..(Q2(IP1 .K)-02( I .K) )-Q2( IM.K) ).V0L 
DI S2( I .3)  —  (Q3( IP2 . K)-3 . • (03( I  PI . K)-03( I , K) )-Q3( IM, K ) ) *VOL 
DIS2( I . 4)  —  (HP2-3 . • (HP  1— HP) -HPM) » VOL 

2  CONTINUE 

DO  15  I  »  1  .  IM1 

D2P  -  AMAX1(EPS(I).EPS(I+1)) 

C22  «  60  .  •  D2P 

Cl  1  -  AMAX 1 (0 . 0 , ( 1 . -C22) ) 

C22  *  C22  •  WW/YACOB( I.K)  •  DT 
C11  -  C11  •  WW/YACOB ( I . K )  •  DT 

C ! ! ! ! ! SWITCH  ON  SECOND-ORDER  AND  SWITCH  OFF  FOURTH-ORDER  DISSIPATION 
C! ! ! ! ! IN  VICINITY  OF  SHOCKS 

DISKI, 1)  -  C11  •  DIS2(  1.1)  +  C22  •  DIS1  ( I  ,  1 ) 

DIS1(I. 2)  -  C11  •  DI S2( I ,2)  +  C22  •  DIS1 (1.2) 

DISKI.  3)  -  C11  *  DIS2(I  .3)  +  C22  «  0IS1(  1.3) 

DISKI. 4)  -  C1I  •  DI$2(  I  ,4)  +  C22  •  DIS1( I .4) 

15  CONTINUE 

DO  16  I  -  2  .  IM1 

DOI(I.K)  -  DISKK 1)  -  DISI(I-KI) 

002 (I.K)  -  DISKI. 2)  -  DISKI-1.2) 

DQ3(  I.K)  -  DIS1(1. 3)  -  DISKI-1.3) 

D04(  I.K)  -  DISK  I  .*)  -  DISKI-1.4) 

16  CONTINUE 
10  CONTINUE 

C  K  DIRECTION 

C ! ! ! ! ! FOURTH-ORDER  DISSIPATION  ONLY 


DO  30  I  -  2  ,  IM1 

WT-  0.5  •  DT  «  WW  /  YACOB( 1,2) 

W3  “  0 . 5  •  DT  •  WW  /  YACOB(I.KMI) 

001(1.2)  -WT*  (01 (1.1)  -  2.  •  01(1.2)  +  01(1.3)) 
1+001(1.2) 

DQ2( 1.2)  -WT.  (02(1.1)  -  2.  »  Q2(1.2)  +  02(1,3)) 

1 +DQ2 (1.2) 

003(1.2)  -WT.  (03(1.1)  -  2.  .  03(1.2)  +  03(1,3)) 

1+DQ3( 1.2) 

DQ4( 1.2)  -WT.  (04(1,1)  -  2.  .  04(1.2)  +  04(1.3)) 

1+DQ4( 1.2) 

WTs  W3 

DO 1(1 .KM1)  -WT.  (01 ( I . KM2)  -  2.  *  Ol(I.KMI)  +  01(1. KMAX)) 
1+DQ1 ( I . KM1 ) 

DQ2 ( I , KM1 )  =WT.  (02(1. KM2)  -  2.  •  Q2(1.KM1)  +  02(1, KMAX)) 
1 +D02 ( I . KM1 ) 
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003  (  I  .  KM  1  )  -WT.  (03  (  I  ,  KM2)  - 
1+003(1 ,KM1 ) 

004(1. KM1)  -WT»  (04( I ,KM2)  - 
1+004(1 .KM1) 

00  35  K  -  3  .  KM2 
WT-  -  DT  •  WW  /  YAC08(  I  ,K) 
001(1. K)  -WT.  (01(1. K+2)  -  4 
1  I  ,K)  -  4  •  01(1. K— 1 )  +01(1 
002 ( I , K )  -WT»(Q2 ( I , K+2)  -  4. 

1  I .K)  -  4.  *02(1 ,K-1 )  +02(1 
003(1. K)  «WT.(Q3(1 .K+2)  -  4. 
II. K)  -  4.  .  03(1, K-1)  +  03(1 
004(1. K)  «WT.(04( I ,K+2)  -  4. 
11. K)  -  4.  .  04(1, K-1)  +  04(1 
35  CONTINUE 
30  CONTINUE 

RETURN 

END 


-  2.  •  03(1. KMI)  +  03(1. KMAX)) 


-  2.  •  04(1, KMI)  +  04(1. KMAX)) 


•  01(1. K+1)  +  6.  •  01! 
K-2))+D01 (I .K) 

•  02(1. K+1)  +  6.  •  02 ( 
K-2) )+DQ2( I .K) 

•  03 ( I ,K+1 )  +  6  •  03( 

K-2 ) )+DQ3( 1 ,K) 

•  04 ( I . K+1 )  +  6.  •  Q4( 
K-2) )+0Q4( I , K ) 


SUBROUTINE  WALLBC 
C0M4ON/SUR  F/PSUR (161) 

COM40N/GRID1A('61  .41  )  ,2(161  .41  ) 

CCM40N/PAR/GAU4A , REYREF . ALFA . ALFA 1  , RED FRE, AM INF, ALFA I 
C0M4ON/DGR ID/DT . I MAX .KMAX . ITEL . ITEU 
COMUtON/GRID/YACOB(161  .41) 

C0M4ON/D  AMP/YRI .  WW 1 

COM40N/MTRIX/XIX(161 . 41  )  ,XI2(  161 . 41 ) , ZETAX( 1 61 .41). 
+2ETA2(161 .41). 

1 X1T( 1 61 .41). ZET  AT (  1 6 1 .41) 

COM40N/FLOW/Q1(161  . 41  )  . 02 (  1 6 1  . 41  )  . Q3(  1 6 1  . 41  )  , 04(  1 61  . 41  ) 
DIMENSION  Cl (4) 

DIMENSION  A(2.2) ,RHS(2) 

C! M ! I  APPLY  BOUNDARY  CONDITIONS  ON  THE  CUT  ANO  THE  AIRFOIL  SURFACE 
DO  9  I-ITEU.IMAX 
II  -  IMAX  +1-1 

01(1,1)  -  .5  •  (01(I.2)-H31(tl  .2)) 

02(1.1)  -5»(Q2(I. 2)+02(I1.2)) 

03(1.1)  -  .5  .  (03(1.2)  +  03(11.2)) 

04(1.1)  -  .5  •  (04(1.2)404(11,2)) 

01(11.1 )— 01 (1.1) 

02(11 . 1 )-Q2( 1.1) 

03(  11,1  )**03(  1,1) 

Q4(I1,1)-Q4(I.1) 

9  CONTINUE 

DO  1  I-  ITEL  .  ITE 
K  -  3 

Cl ( 1 )  -  XIT(I.K) 

Cl  (2)  -  XIX(I.K) 

Cl  (3)  »  XIZ(I.K) 

UC0N3  -  (Q2(T.K)«C1(2)+03(I.K)»C1(3)) 

1/01 ( I . K ) 

K-2 

Cl ( 1 )  -  XIT(I.K) 

Cl (2)  -  XIX(I.K) 

Cl  (3)  -  XIZ(I.K) 

UC0N2  -  (02( I ,K)*C1 (2)+03( I , K ) *C 1(3)) 

1/01(1. K) 

RHS( 1 )  -  2.  •  UC0N2  -  UC0N3  -  XIT(I.I) 

C  FOR  VISCOUS  FLOWS  SET  UCON  TO  ZERO  ALSO 
I F ( REYREF . GT .0. )  RHS(1)  -  -  XIT(I.I) 

A(  1  . 1  )  -  XIX(1 . 1) 

A(1.2)  -  XIZ(I.I) 

A(2,1)  -  ZETAX(I.I) 

A(2 , 2)  -  ZETAZ(I.I) 


ooo  oo  o  o  o 
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RHS ( 2 )  -  -  ZETAT (1,1) 

TEMPI  -  A( 1  .  1 ) 

TEVP2  -  A(1 ,2) 

TEMP3  -  A(2, 1 ) 

TEMP4  -  A(2. 2) 

DEN  -  t.  /(TEMPI  .  TEMP4  -  TEMP2  •  TEMP3) 

A(1 ,1)  -  A(2 .2)  •  DEN 

A(1 .2)  -  -  TEMP2  •  DEN 

A(2 . 1 )  -  -  TEMP3  •  DEN 

A(2 , 2)  -  TE3k^1  •  DEN 

01(1.1)  -  2.  •  01(1.2)  -  01(1,3) 

02(1.1)  -  Q1(I.1)«(A(1.1 ) »RHS( 1 )+A( 1 , 2) «RHS( 2) ) 
03(1,1)  -  01(1.1)  • ( A(2 . 1 ) «RHS( 1 )+A( 2 . 2) »RHS(2 ) ) 

1  CONTINUE 

DO  10  1-ITEL  , ITEU 
U2=02( 1 . 2 )/01 ( I .2) 

W2“03( I . 2)/01 (1,2) 

P2»(GAM*A-1  .  )«(Q4(  I  ,2)-0.5.Q1(  1 , 2) • (U2«U2+W2*W2) ) 
U3-02(I,3)/01(I.3) 

W3-Q3(I.3)/01(I.3) 

P3-(GAAMA-1  •  )»(Q4(I  ,3)-0.5*Q1(I  .3)*(U3.U3+W3«W3)! 
Pl«(4. «P2-P3)/3 . 

PSUR(  1  )«(GAjy*UA»Pl  — 1  .  )/(  .  7«AMINP *»2) 

Ul«02( 1 . 1 )/Q1  (1.1) 

#1-Q3( I . 1 )/Q1 (1.1) 

10  04(1.1 )-P1/(GAM4A-1 . )+0.5»O1 ( I . 1 ) • (U1 «U1+W1 «W1 ) 
RETURN 
END 


SUBROUTINE  STRESS(ITN.DALFA) 

COM4ON/FLOW/01 (161 .41 ). 02(1 61 ,41 ) ,03( ’61 .41) .04(161 .41) 
C0M4ON/DGRID/DT , I MAX .KMAX , I  TEL . ITEU 
C0M40N/GRID1/X(161 ,41).Z(161 .41) 

COMMON/P  AR/GAM4A .  REYREF , ALFA  . ALFA1  .REDFRE .  AMINF  ,  ALFAI 
C0M40N/PERTR/DQ1 (161 . 41 ) .002 ( 1 61 . 41 ) .DQ3(161 .41 ) .D04(161 .41) 
C0M40N/MUTUR/CMU(161 .41) 

DIMENSION  AA( 1 61 .41) 

1.RH 1(161). RH2 (161) . RH3( 161 ) ,RM4(161) 

COMMON/LOG IC/RSTRT .PITCH .PLUNGE 
LOGICAL  RSTRT. PITCH. PLUNGE 
U(I.J)  -  02(1. J)  /  01(1. J) 

V(I.J)  -  03(1. J)  /  01(1. J) 

THIS  SUBROUTINE  ADDS  VISCOUS  TERMS  TO  THE  RIGHT  HAND  SIDE 
GOGM  -  GAMM  /  (GANMA  -  1  ) 

I F( I TN . GT . 10. OR . (RSTRT) )  CALL  EDDY(DALFA) 

COMPUTE  U  AND  V 
KMAXM1  -  KMAX  -  1 
IMAXM1  -  IMAX  -  1 
PR  -  1  . 

DO  10  K  *  1  .  KMAX 
DO  10  I  >  1  ,  IMAX 

E  -  04(1. K)  /  01(1, K)  -  0.5  •  (U( I , K ) • • 2+V ( I , K ) • • 2 ) 

10  AA( I . K )  -  GOGM  •  E 

COMPUTE  TXX.TXY  AND  VISCOUS  DISSIPATION  AT  I  -  1  /  2 

DO  30  K  -  2  .  KMAXM1 
DO  20  I  -  2  .  IMAX 
UXI  »  U(I .K)  -  U( 1-1 ,K) 

VXI  -  V(I.K)  -  V(l-I.K) 

AX I  -  AA( I .K)  -  AA(I-I.K) 

UZET-  ,25»(U( I ,  K+1 )— U( I .K-1 )+U( 1-1 ,K+1 )-U( 1-1 ,K-1 ) ) 

VZET«  . 25» (V( I . K+1 )-V( I .K-1 )+V( 1-1 .K+1 )-V( 1-1 .K-1 )) 

AZET-  . 25* ( AA( I .K+1 )-AA(I .K-1 )+AA(I-1 ,K+1)-AA(I-1 ,K-1)) 

XXI  -  X ( I , K)  -  X(I-I.K) 
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ZXI  -  Z( I . K)  -  Z( 1-1 .K) 

XZET-  .25  •  (X( 1 ,K+1 ) — X ( I ,K— 1 )+X( 1-1 ,  K+1 )-X( 1-1 ,K— 1 ) ) 
ZZET-  .25  •  (Z( I .K+1 )-Z( I . K— 1 )+Z( 1-1 ,  K+1 )-Z( 1-1 ,K— 1 ) ) 

YAC  -  XXI  •  ZZET  -  ZXI  •  XZET 
VAC  -  1 .  /  YAC 
XIX  -  ZZET  •  YAC 
ZETAX-  -  ZXI  .  YAC 
XIZ  -  -XZET  •  YAC 
ZETAZ-  XXI  •  YAC 

CNM  -  .5  •  (CMU(I.K)  +  CMU(I-I.K)) 

UX  -  UXI  •  XIX  +  UZET  .  ZETAX 
VX  -  VXI  •  XIX  +  VZET  •  ZETAX 
AX  -  AX  I  •  XIX  +  AZET  .  ZETAX 

UZ  «  UXI  .  XIZ  +  UZET  .  ZETAZ 

VZ  -  VXI  «  XIZ  +  VZET  .  ZETAZ 

AZ  -  AX  I  •  XIZ  +  AZET  .  ZETAZ 

TXX  -  -(-4.  .  UX  +  2.  •  VZ)  •  CNM  /  3. 

TXY  -  CNM  •  (UZ  +  VX) 

TYY  «  -CNM  /  3.  •  (-4.  •  VZ  +  2.  •  UX) 

R4  -  ( ( U( I , K  )+U( 1—1 ,K))«TXX+(V(I .K )+V( 1—1 , K ) ) »TXY) *0 . 5 
1  +  CNM  /  PR/(GAN*4A  -  1 . )  •  AX 

S4  -  ( (U( I ,K)+U( 1-1 ,K) ).TXY+(V( I .K)+V( 1—1 .  K )  ) «TYY) «0 . 5 
1  +  CNM  /  PR  /  (GAVI4A  -  1  .  )  •  AZ 

DEBUG 

TURN  OFF  ENRGY  DISSIPATION  AND  DIFFUSION 
R4  -  0. 

S4  »  0. 

RH 1(1)  -  0. 

RH2( I )  -  (XIX  •  TXX  +  XIZ  •  TXY)  /  YAC 

RH3( I )  -  (XIX  .  TXY  +  XIZ  •  TYY)  /  YAC 

20  RH4( I )  -  (XIX  .  R4  +  XIZ  •  S4)  /  YAC 

DO  30  I  -  2  .  IMAXM1 

DQi(I.K)  -  DQI(I.K)  +  RH1 ( 1+1 )  -  RH1(I) 

OQ2(I. K)  -  D02 ( I . K )  +  RH2(I+1)  -  RH2( I ) 

DQ3(  I  .  K )  -  003(1. K)  +  RH3(I-H)  -  RH3(I) 

D04 ( I , K )  -  D04 ( 1 . K )  +  RH4( 1+1 )  -  RH4(I) 

30  CONTINUE 
C  IN  THE  Z  DIRECTION 

DO  70  1  -  2  .  IMAXM1 

DO  60  K  *  2  .  KMAX 

UXI  -  .25  •  (U(I-t-I.K)— U(I-1.K)+U(I  +  1.K-1)-U(I-1.K-1)) 

VXI  -  .25  •  ( V( 1  +  1 , K )— V ( I  — 1 , K )+V ( 1  +  1 , K— 1 )— V FI  —  1 ,X— 1  ) ) 

AX  I  -  .25  •  ( AA( 1  + 1  , K )— AA( I  —  1 , K )+AA( 1  +  1 , K— 1 )— AA( I— 1 . K— 1 ) ) 
XXI  -  .25  •  ( X ( I +1 , K )— X ( I— 1 . K )+X ( I +1 . K— 1 )— X (1—1 ,K— 1 ) ) 

ZXI  -  .25  •  (Z( 1  +  1 ,K)-Z( 1-1 ,K)+Z(I  +  1 ,K-1 ) — Z (I  —  1 ,K— 1  )) 

UZET  -  U(I.K)  -  U(I.K-I) 

VZET  -  V( I .K)  -  V( I , K-1 ) 

AZET  -  AA(I.K)  -  AA(I.K-I) 

XZET  -  X(I . K)  -  X(I . K— 1 ) 

ZZET  -  Z(I .K)  -  Z(I . K-1 ) 

YAC  -  XXI  •  ZZET  -  ZXI  •  XZET 

YAC  -  1  /  YAC 

XIX  -  ZZET  «  YAC 

ZETAX-  -  ZXI  »  YAC 

XIZ  -  -XZET  •  YAC 

ZETAZ-  XXI  •  YAC 

CNM  -  .5  •  (CMU( I .K)  +  CMU( I ,K-1 )) 

UX  -  UXI  •  XIX  +  UZET  .  ZETAX 
VX  -  VXI  •  XIX  +  VZET  .  ZETAX 
AX  -  AX  I  •  XIX  +  AZET  •  ZETAX 

UZ  -  UXI  •  XIZ  +  UZET  .  ZETAZ 

VZ  -  VXI  •  XIZ  +  VZET  .  ZETAZ 

AZ  -  AX  I  •  XIZ  +  AZET  .  ZETAZ 

TXX  -  -(-4  .  UX  +  2.  •  VZ)  •  CNM  /  3. 

TXY  -  CNM  «  (UZ  +  VX) 

TYY  -  -CNM  /  3.  •  (-4.  .  VZ  +  2.  •  UX) 

R4  -  ( (U( I  .K)+U(l .K-1 )).TXX+(V( I .K)+V(I . K-1 ) ) .TXY) »0 . 5 
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1  +  CNM  /  PR/(GA)*4A  -  1 . )  •  AX 

S4  -  (  (U(  I , K  )+U( I . K— 1 ) ) »TXY+( V( I , K  )+V ( I , K— 1 ) ) • TYY )  »0 . 5 

1  +  CNM  /  PR  /  (GAMMA  -  1 . )  •  AZ 

R4  -  0. 

S4  -  0. 

RH 1 ( K )  -  0. 

RH2 (K )  -  (ZETAX  .  TXX  +  ZETAZ  •  TXY)  /  YAC 

RH3(  K )  -  (ZETAX  .  TXY  +  ZETAZ  •  TYY)  /  YAC 

60  RH4(K )  -  (ZETAX  .  R4  +  ZETAZ  •  S4)  /  YAC 

DO  70  K  -  2  .  KMAXM1 

DQI(I.K)  -  001(1. K)  +  RHI(K-M)  -  RH1(K) 

D02 ( 1 , K )  -  D02(I.K)  +  RH2 ( K+1 )  -  RH2(K) 

DQ3(1.K)  -  D03( I ,K)  +  RH3(K+1)  -  RH3(K) 

D04 ( I . K )  -  DQ4( I ,K)  +  RH4(K+1)  -  RH4(K) 

70  CONTINUE 

RETURN 

END 


C 


C 


SUBROUT  I NE  LOAD (CPS . CL . CD . CM , ALFAS ) 
C0M40N/GRID1/X(161 . 41 ) . Y ( 1 6 1 .41) 

C0M40N/DGR I D/DT , I MAX . KMAX . I T  E  L . ITEU 
DIMENSION  CPS( 161 ) 

THIS  SUBROUTINE  COMPUTES  THE  INVISC1D  CONTRIBUTIONS 
TO  LOADS  ON  THE  AIRFOIL  SURFACE 

IMAXM1  -  IMAX  -  1 
CL  -  0. 

CD  -  0. 

CM  -  0. 

DO  400  I  -  ITEL  .  ITEU  -  1 
XL  -  .5  •  (X(I,1)+X(I+1,1)) 

YL  -  5  •  ( Y ( I . 1 ) +Y (1+1,1)) 

DX  -  X(l+1 .1)  -  X(I.I) 

DY  .  Y( 1+1 , 1 )  -  Y( I . 1 ) 

CPA  -  CPS(I+1)  •  .5  +  CPS(I)  •  .5 
DCL  «  CPA  .  (-0X) 

DCD  «  CPA  .  DY 
CL  -  CL  +  DCL 
CD  -  CD  +  DCD 

400  CM  -  CM  +  DCD  •  YL  -  DCL  •  XL 

DCL  -  CL  •  COS(ALFAS)  -  CD  »  SIN(ALFAS) 

CD  -  CL  •  SIN(ALFAS)  +  CD  •  COS(ALFAS) 

CL  -  DCL 

RETURN 

END 


SUBROUTINE  WRAP(II , JJ.XSING. YSING.XP.YP.S0.A0.Y0) 

DIMENSION  S0( 1 61 .4) , Y0(41 . 4) . A0( 1 61 . 4) , XP( 1 ) , YP( 1 ) 

THIS  SUBROUTINE  UNWRAPS  THE  AIRFOIL  AND  STORES  THE  UNWRAPPED 
SURFACE  IN  ARRAYS  A0  AND  S0 .  IT  ALSO  DETERMINES  THE  STRETCHING 
IN  THE  ETA  DIRECTION. 

IMID  -  (II  +  1)  /  2 
DY  -  8  /  (JJ  -  2) 

DO  1  J  -  2  .  JJ 
Y  -  FLOAT ( J— 2 )  •  DY 
1  Y0(J,1)  -  1.25  •  Y  /  (1.  -  Y  .  Y) 

Y0( 1  ,  1 )  -  -  Y0(3, 1) 

PI  -  4.  •  ATAN  (  1 . ) 


oooo  ooo  o  ooo  ooo 


PI 

XSING 

YSING 


-  1 


II 

-  XSING 

-  YSING 

+  ATAN2( (U»Y1 1— V»X1 1  ' 
+  Y11..2) 


,  (U«X1 1+V»Y1 1 ) ) 


ANGL  -  PI 
U  -  XP(1) 

V  «  YP(1) 

U  -  1  . 

V  -  0. 

I  I  Ml  -  II 
DO  2  I  -  1  . 

XII  -  XP( I ) 

rn  .  yp(  I ) 

ANGL  -  ANGL 
R  -  SORT (XI 1 »»2 
U  -  XI  1 

V  -Y11 
R  -  SORT(R) 

A0(1.1)  -  R  .  COS ( . 5 

2  S0( 1.1)  -  R  •  SIN( .5 
!  !  ! !  ! I r  OUTPUT  OF  UNWRAPPED 
WRITE  (6.1000) 

WRITE  (6.1000)  (I,A0(I.1),S0(I,f).I  - 
RETURN 

FORMAT (IX.- UNWRAPPED  COORDINATES  IN  THE  TRANSFORMED  PLANE’) 
FORMAT (15  .  2F20.8) 

END 


ANGL) 

ANGL) 

COORDINATES 


IS  DESIRED 


f 


II) 


1000 

2000 


SUBROUT  IN:.  TAB  I  NT(  XP .  YP ,  XS I NG .  YS I NG .  N) 

DIMENSION  XP(161).YP(l6t),S0(161),A0(161) 

!  ! ! ! ! SMOOTH  THE  AIRFOIL  SURFACE  BY  FINDING  ADDITIONAL  POINTS 
U  -  XP(1)  -  XSING 

V  -  YP( 1 )  -  YSING 
U  -  1  . 

V  -  0. 

ANGL  -  8  •  ATAN( 1 . ) 

DO  1  I  -  1 . N 

XII  -  XP( I)  -  XSING 

Y1 1  -  YP( I)  -  YSING 

ANGL  -  ANGL  +  ATAN2( (U»Y1 1-V.X1 1 ) . (U»X1 1+V»Y1 1 ) ) 

R  *  SORT(X7 I*»2  +  Y1>  •«  2) 

U  -  XII 


•5) 

•5) 


V  -  Y11 
R  -  SORT (R) 

A0(I)  -  R  .  COS ( ANGL 
S0(I)  -  R  •  S I N( ANGL 
DX  -(A0(N)-A0(1 ))/96. 

A0ST  -  A0( 1 ) 

DO  3  I  ■  1  .97 

XX  «  FLOAT ( 1-1 )  »  OX  +  A0ST 

CALL  TAINT (A0. S0. XX. YY. N. 0, NER.MON) 

XP(I)  -  XX  •  XX  -  YY  .  YY  +  XSING 

YP( I )  -  2 .  •  XX  «  YY  +  YSING 

RETURN 

END 


SUBROUTINE  TAINT(XTAB . FTAB , X . FX , N , K , NER .MOf  ) 
DIMENSION  XTAB(1).FTAB(1).T(10),C(10) 


NASA  -  AMES  SUBROUTINE  FOR  POLYNOMIAL  INTERPOLATION 
OF  A  TABULATED  FUNCTION 


IF(N-K)  1  ,  1 

1  NER  -  2 
RETURN 

2  I F ( K— 9)  3.3.1 

3  IF(MON)  4,4.5 


3 


IF(M0N-2)  6.7.4 
J  -  0 

NM1  -  N  -  1 

DO  8  I  -  1  .  NM1 

IF(XTAB(I)  -  XTAB( 1+1 ) )  9.11.10 

NER  -  3 

RETURN 

J  -  J-1  . 

GO  TO  8 
J  «  J  +  1 
CONTINUE 
MON  -  1 

IF(J)  12  .  6  .  6 
MON  -  2 

DO  13  I  -  1  .  N 

I F(X  -  XTAB(I))  14.14.13 

J  *  I 

GO  TO  18 

CONTINUE 

GO  TO  15 

DO  16  I  -  1  ,  N 

I F ( X-XTAB( I ) )  16.17,17 

J  -  I 

GO  TO  18 

CONTINUE 

J  -  N 

J  -  J  -  (K+1 )  /  2 

I F( J )  19,19.20 

J  -  1 

M  -  J  +  K 

I F(M  -  N)  21 .21 .22 

J  =  J  —  1 

GO  TO  20 

KP 1  -  K  +  1 

JSAVE  -  J 

DO  23  L  -  1 .  KP1 

C(L)  -  X  -  XTAB(J) 

T(L)  -  FTAB(J) 

J  -  J+1 
DO  24  J  «  1.K 
I  •»  J+1 

T ( I )  -  (C(J).T(I)-C(I).T(J))/(C(J)-C(I)) 
I  -  1  +  1 

I F ( I-KP1 )  25.25.24 

CONTINUE 

FX  -  T (KP1 ) 

NER  -  1 

RETURN 

END 


1  ,  KP1 
XTAB(J) 


Z2  ••  2  -  21  ••  2 

X2  -  XI 
22  -  Z1 

X3  ••  2  —  XI  ••  2 
Z3  ••  2  —  Z1  •  •  2 

X3  -  XI 
Z3  -  21 

(D7  •  (  D1  +  D2)  - 


B  -  (D7  •  (  D1  +  D2)  -  D3* (D5+D6) )/(2  »(D7»D4— D3»D8) ) 
I F ( ABS(D3) . LT . ABS(D7 ) )  GO  TO  10 
A  -  (D1  +  D2  -  2.  •  B  •  04)  /  (2.  •  03) 


GO  TO  20 

10  A  -  (05  +  06  -  2.  •  B  •  08)  /  (  2.  •  07) 
20  CONTINUE 

R  -  SORT ( (X2-A) • •  2  +  (Z2-B)»«2) 

XLE  -  X ( NU) 

YLE  -  Z(NU) 

CHD  -  X( 1 )  -  X(NU) 

A2  -  ( Z ( 2 )— Z ( 1 ) )  /(X(2)  -  X( 1 ) ) 

A1  -  ( Z (N)— Z(N-1 ) )/ ( X(N)— X(N— 1 ) ) 

TES  -  .5  •  ( A 1  +  A2) 

TEA  -  A2  -  A1 
TEA  -  TEA  «  57.29578 
XSING  -  (A+XLE)  /2 . 

YSING  -  (B+YLE)  /  2. 

RETURN 
ENO 


SUBROUTINE  AIRFOL( 1 1 . JJ . TT . IE) 

C0U*)N/GRID1/X(  161 .41  )  .  Z  ( 1  61  ,41) 

CCA4AON  /YS  YM/ 1 SYM 

DIMENSION  S0(161 .4) .A0(161 .4) .Y0(41 . 4) . XP( 1 61 ) . YP( 1 61 ) . 
1E(161),F(161).B0(49) 

DATA  (B0(1). 1-1 ,32)/1 . .1 .0414. 1 .0836.1  1270.1. 1715.1.2175,1.2651 . 
11 .3145.1.3659.1 .4199.1 .4755.1.5349.1.5973.1-6636.1.7342.1.8099. 

21 .8914. 1 .9799.2.0764.2. 1829.2.3012.2.4341 .2.5653.2.7597.2.9646. 
33.2106.3.5141 .3. 90 19. 4. 42 19. 5. 1687.6.3632.8.6809/ 

! ! ! ! COMPUTE  THE  COMPUTATIONAL  GRID  POINTS  BASED  ON  INPUT  AIRFOIL  SHAPE 
DO  8  I  -  1  .32 
8  A0 (  I  .  1 )  -  B0(I) 

READ  (5,1) 

1  FORMAT (IX) 

READ  (5.2)  FNU.FNL.ZSYM 

2  FORMAT ( 3F1 0 . 0) 

I SYM  -  0 

I F ( ZSYM . NE ■ 0 . )  I SYM  -  1 
II  -  157 
JJ  «  41 

IT  -  31 

IE  -  127 
I  I  PI  -  II  +  1 
I  IM1  -  I  I  -  1 
I  I JJ  -  I  I  •  JJ 
I IJJ2  -  II  •  (  JJ-2) 

ILE  -  (IT  +  IE  )  /  2 

ISTP  -  0 

NN  -  5 

NRF  -  0 

NOTAPE  -  1 

PI  -  4.  «  ATAN( 1 . ) 

NU  »  FNU 
NL  -  FNL 
N  -  NU  +  NL  -  1 
READ ( 5 . 1 ) 


TOWS 


I 


o  o  o  o  o  o  o 


mwuuviwuwnu!  wtowottotowttto  *f  ilwjuwjvjvw.wata' w 


I 


I 


! 


I 


READ  (5,333)  (XP( I ) . YP( I ) . J  -  NL  ,  N) 

333  FORMAT ( 2F 1 0 . 0) 

9994  FORMAT (F20. 8) 

L  -  N  +  1 

I F ( ZSYM  GT.  0. )  GO  TO  9995 
L  -  NL  +  1 
READ(5 , 1 ) 

READ  (5,333)  ( XP( L-l ) . YP( L-l ) . 1-1 , NL) 

GO  TO  9996 

9995  K1  -  L 

DO  16  I  -  NL  .  N 
K  -  K1  -  I 
XP(K)  -  XP(  I  ) 
vp(K)  -  -  YP( I ) 

16  CONTINUE 

9996  SCALE  -  1-  /(XP(I)-XP(NL)) 

XX  «=  XP(NL) 

YY  =  YP(NL) 

DO  9997  I  =  1  .  N 
XP(I)  «  XP(I)  •  SCALE 

9997  YP( I )  -  YP( I )  •  SCALE 

CALL  S I NG ( NU . N . XP , YP . XLE , ZLE , TEA , TES , XS I NG , YS I NG , CHD ) 

CALL  TAB  I  NT ( XP , YP , XS I NG , YS I NG , N) 

NBODY  -  IE  +  1  -  IT 
DO  6791  1-1  .  NBODY 
L  «  I  -  1 
E( IT+L)  =  XP(I) 

6791  F( IT+L)  -  YP ( I ) 

IEP1  -  IE  +  1 
SLOPT  =  TES  •  ,75 
DO  438  I  -  IEP1  .  II 
II  -  I  +1  -  IE 
E ( I )  -  A0 (11,1) 

DXI  -  1 .  /  48. 

D  -  4.  /  3.  •  ( E( I )  -  .25) 

F ( I )  —  F(IE)  +  SLOPT  .  ALOG(D)  /  D 
L  -  IIP1  -  I 
E(L)  -  E( I  ) 

438  F(L)  -  F ( I T )  +  SLOPT  .  ALOG(D)/D 

C  INRITE  (6,439) 

439  FORMAT ( 2X , 3H  I , 19X , 1HX, 19X . 1HY) 

C  WRITE  (6,37)  (I.E(D.F(I).I  -  1  .  II) 

CALL  WRAP ( I  I , J  J , XS I NG . YS I NG , E . F . S0 , A0 , Y0 ) 

C! ! ! ! IMAP  GRID  BACK  TO  PHYSICAL  PLANE  AND  SHIFT  TO  QUARTER  CHORD 
DO  10  J  -  2  .  JJ 
DO  10  I  -  1  ,  II 

X(I.J-I)  ■=  A0(I,1)..2  -  (S0( I . 1 )+Y0( J , 1 ) ) »»2 
1  -  0  25 

10  Z(I.J-I)  -  2.  •  A0(I,1)  .  (S0( I . 1 )+Y0( J , 1 ) ) 

JJ  -  JJ  -  1 
RETURN 

37  FORMAT ( I5.2F20.8) 

END 


SUBROUTINE  CLUSTR(DS) 

COMMON/GR I D 1 /X (161  ,  41 ) .2( 161 ,41 ) 
CONMON/DGRID/DT, I MAX , KMAX , I  TEL. ITEU 
DIMENSION  S(41),XP(41).YP(41),R(41) 


THIS  SUBROUTINE  CLUSTERS  A  GIVEN  X.Z  GRID  SUCH  THAT  THE  FIRST  POINT  IS  AT 
THE  USER-SPECIFIED  DISTANCE  DNMIN 
!  !  !  !  (COMPUTE  THE  OLD  STRETCHING 
DO  100  I  -  1  .  IMAX 
S(1)  -  0 
XP(1 )  -  X(  I  ,  1  ) 


yp(i) 

DO  10  K  -  2  ,  KMAX 
XP(K)  -  X( I ,  K) 
yp(K)  -  Z( I ,K) 

10  S(K)  -  SORT((XP(K)-XP(K-1)).*2+(YP(K)-YP(K-1))..2) 
1+S(K-1) 

SUMDX  -  S ( KMAX ) 

CALL  STRTCH( SUMDX, OS. H .KMAX , FACTOR) 

C  WRITE  (6,200)  I. FACTOR 

R(1)  -  0. 

DR  -  DS 

00  20  K  -  2  ,  KMAX 
R(K)  -  R(K-1 )  +  DR 
DR  -  DR  •  FACTOR 
20  CONTINUE 

RLAST  -  1 .  /  R(KMAX) 

DO  30  K  «  2  .  KMAX 
R 1  -  R(K)  •  RLAST  •  S(KMAX) 

CM  ! ! ! REDISTRIBUTE  THE  CONSTANT-ETA  LINES 

CALL  TAI NT ( S . XP . R 1 ,XP1 . KMAX . 3 . NER , MON) 

X( I ,K)  -  XP1 

CALL  TAINT(S. YP.R1 . Y PI , KMAX , 3 , NER . MON) 

Z(I  .K)  -  YP1 
30  CONTINUE 
100  CONTINUE 
C  WRITE  (6.115) 

DO  1 1 0  I  -  1  .  I  MAX 
DX  -  X( I .2)  -  X(I.I) 

OY  =.  Z(  I  .  2)  -  Z(I.I) 

DN  -  SQRT(DX«DX+DY.DY) 

C  WRI TE( 6 . 120)  I  .  DX  .  DY  .  DN 

110  CONTINUE 
RETURN 

115  FORMAT (5X.6HNORMAL, IX. 8HD I  STANCE. 3H  AT, AH  THE.5H  WALL./ 
1 ,5H  I . 8X . 2HDX . 8X . 2HDY ,8X . 2H0N ,//) 

120  FORMAT ( I  5 . 3F1 0  5) 

200  FORMAT (I5.F10.5) 


END 


SUBROUTINE  STRTCH( SUMDX .0X1 .FI ,N1 .R) 

THIS  SUBROUTINE  DETERMINES  A  GEOMETRIC 
PROGRESSION  OF  GRID  SPACING  BETWEEN  1  AND  N1  SUCH  THAT 
SUM8DX)  EOUALS  SUMDX.  THE  RATIO  BETWEEN  SUCCESSIVE 
SPACINGS  IS  R. 

N  -  N1  -  1 
R  -  1  .5 
El  -  1 . E-04 
E2  -  1 • E-04 
DO  1 0  L  -  1 ,  50 

F-  (R-1 )  .  SUMDX  -  DX1.(R..N-1) 

FP  =  SUMDX  -  DX1  •  FLOAT (N)  •  R  ••  (N-1) 

RITER  -  R  -  F/  FP 

IF(1 .E-02.LT. RITER. AND. RITER. LT.1 . )  RITER  -  1. 

I F(  1  .  .LT. RITER. AND. RITER.LT. 100.)  RITER-. 01 
I F( ABS(R-RI TER) . LT .  R»E 1 )  GO  TO  1 
R  -  RITER 
10  CONTINUE 
R  -  1 .0001 

DX1  -  DZT0T/FL0AT(N1-1) 

RETURN 
1  R-  RITER 
RETURN 


c 


SUBROUTINE  EDDY ( OALFA) 

COMAON/FLOW/Q1 ( 1 6 1 . 4 1 ) . Q2 ( 1 6 1 . 41 ) . 03 ( 1 6 1 , 41 ) . 04 ( 1 6 1 . 4 1 ) 
COMMON/MU TUR/CMU( 161 ,41  ) 

C0M40N/SK I NC F/C F ( 1 6 1 ) 

COMMON/DGR I D/DT , IMAX.KMAX.  I  TEL. ITEU 

COMUON /PAR/GAVNA  ,  REYREF  ,  ALFA ,  ALFA  1  .  REDFRE  ,  AMI  NF  .  ALFA  I 
COMMON/GR I D 1 /X ( 1 6 1 ,41),Z(161,41) 

Dimension  tin(4i ) ,tout(4i ) , y(ai ) 

INITIALIZE  VISCOSITY  EVERYWHERE 
FACT1  -  DT  •  AMINF  /  REYREF 
CMUMAX  -  100.  •  FACT  1  /  DT 
DO  1  K  -  1  .  KMAX 
DO  1  1-1  .  I  MAX 
1  CMU( I ,K)  »  FACT  1 

THIS  SUBROUTINE  COMPUTES  THE  EDDY  VISCOSITY  BASED  ON  THE 

baldwin-lomax  two  layer  model 

ro  100  I  -  2  ,  IMAX  -  1 
UDIF  -  0. 

FMAX  -  0. IE-06 
YMAX  -  . 1 E-06 
FYMAX  -  0. 

Y ( 1 )  -  0. 

UWALL  -  0- 

I F ( I GT . ITEU.OR. I .LE. ITEL)UWALL  -  SORT (02(1,1 )»«2+Q3( I . 1  ).«2)/ 
101(1,1) 


COMPUTE  TAU  AT  THE  WALL 
UET  -  1  • f 02 ( I ,2)/Ql (1.2) 
VET  -  1  «(03( I ,2)/Ql (I .2) 
XXI  -  X( 1  +  1  .  1  )  -  X( I -1 . 1  ) 
ZXI  -  Z( 1+1,1)  -  Z(I-1 . 1) 
XET  -  4.  •  X( I ,2)  -  3.  • 
ZET  -  4  •  Z( I .2)  -  3.  • 


VET  -  1 . 
XXI  -  X( 
ZXI  -  Z( 
XET  -  4 
ZET  -  4. 
XXI  -  .6 
ZXI  -  .5 
XET  -  ,f 
2ET  - 
YAC  -  1 . 
OMEGA  » 
TWAUL  - 
CF( I )  - 
FACT  -  S 
DO  10  K 


02 ( I , 1 )/Q1 (1,1)) 

03(1,1 ) /O 1 ( I , 1 )) 


•  X(I,1)  -  X( I .3) 
.  Z ( I . 1 )  -  Z( I , 3) 


.5  •  XXI 
.5  •  ZXI 
.5  •  XET 
.5  •  ZET 

1 .  /  (XXI  •  ZET  -  ZXI  •  XET) 

-  (UET  .  XXI  -  VET  •  ZXI  )  *  YAC 

-  AMINF  .  OMEGA  /  REYREF 

-  2  •  TWALL  /  (AMINF. .2) 

■  SQRT(Q1(I.1)  .  ABS(TWALL))«REYREF/(26. .AMINF) 


UXI  -  (02(1+1 . K)/01 ( 1+1 ,K)  -  02(1-1 ,K)/Q1(J-1,K)) 

VXI  =  (03(1+1 . K)/Q1 ( 1+1 . K )-03( 1-1 ,K)/01(I-1 ,K)) 

UET  =  (02( I , K+1 )/01 ( 1 , K+1 )-Q2( I . K-1 )/Q1 ( I , K-1 ) ) 

VET  -  (03(1. K+1)/Q1(I.K+1)-03(I.K-1)/Q1(I, K-1)) 

XXI  -  X( 1+1 ,K)  -  X( 1—1 ,K) 

ZXI  -  Z( 1+1 ,K)  -  Z( I— 1 , K ) 

XET  =  X( I , K+1 )  -  X ( I , K— 1 ) 

ZET  -  Z( I , K+1 )  -  Z(I , K— 1 ) 

YAC  -  1 .  /  (XXI  .  ZET  -  ZXI  *  XET) 

OMEGA  -  ABS( UET .XXI+VET »ZX I-UX I .XET-VXI .ZET )  .  YAC 
UDIF  -  SOPT (02(1 , K ) • «2+Q3( I.K)««2)/Q1(I,K)  -  UWALL 
I F ( A8S(UD I F ) . GT . UDI FMAX )  UDI FMAX  -  ABS(UDIF) 

Y(K)  -  SORT ( ( X ( I , K )— X( I ,K— 1))..2+(Z(i ,  K  )  — Z  (  I , K— 1 ) ) * .2 )+Y ( K— 1 ) 

F  =  Y(K)  •  OMEGA 

IF(( y(k).FACT) .GT.20. )  GO  TO  31 

I F ( I  GT. ITEL.AND. I .IT. ITEU)  F  -  F  .  (1.  -  EXP(-Y(K ).FACT ) ) 

31  CONTINUE 

MODIFIED  TUR8ULENCE  MODEL  APPLY  FOR  SPECIFIED  RANGE  OF  ANGLES  WHERE 
FY  IS  USED  T  0  F I  NO  THE  SECOND  PEAK  VALUE  OF  F  FUNCTION 


I F(ALFA . LT . ALFAI . AND . DALFA . GE . 0 . )  THEN 
FY  -  F  «  Y(K) 

IF(FY.GT.FYMAX)  THEN 
FYMAX  -  FY 
FMAX  -  F 
YMAX  -  Y(K) 

END  IF 
END  IF 

I F( ALFA . GE . ALFAI . OR . DALFA. LT . 0 . )  THEN 
IF(F.GT.FMAX)  THEN 
FMAX  -  F 
YMAX  -  Y(K) 

END  IF 
END  IF 

FCT  -  Y(K)  •  FACT 
I F(FCT . GT . 20 . )  FCT  -  20. 

FCT  -  ABS(FCT) 

EL  -  .4  .  Y(K)  •  (1 ■  -  EXP(-FCT)) 


TIN(K)  «  01 ( I .K)  •  EL  •  EL  •  OMEGA 


T I N(K  - 
10  CONT 
KSWT''  ■ 

FWAKE  -"MAX  •  FMAX 

FI  -  •  :5  •  YMAX  •  UDIF  »»2  /  FMAX 

IF(F1  _ 7. FWAKE)  FWAKE  -  FI 

DO  20  X  -  2  ,  KMAX  -  1 

FKLEB  -  0. 

IF(ABS(Y(K)/YMAX).LT. 1 .E+04)  THEN 

FKLEB  -  1.  /  (1.  +  5.5  •  (0.3  •  Y (K )/YMAX )  ••  6) 

END  IF 

TOUT (K)  -  .0168  •  1.6  •  01(1, K)  •  FWAKE  •  FKLEB 
TOUT (K )  -  A8S( TOUT ( K ) ) 

I F (KSWTCH  NE  0)  GO  TO  20 
IF(TIN(K)  GT. TOUT (K ) )  KSWTCH  -  K  -  1 
20  CONTINUE 

CM!!! TOTAL  VISCOSITY  IS  SUM  OF  LAMINAR  AND  INNER/OUTER  LAYER  AS  APPROPRIATE 
DO  30  K  -  2  .  KMAX  -  1 
I  F(K.LE. KSWTCH)  THEN 

CMU(I.K)  -  DT  .  (AMINF/REYREF  +  ABS(TIN(K))) 


A8S(TIN(K) ) 


CMU(I.K)  -  DT  .  (AMINF  /  REYREF  +  ABS(TOUT(K) ) ) 
END  IF 


30  CONTINUE 
100  CONTINUE 
RETURN 
END 


SUBROUTINE  RESI (OMEGA) 

CCAW0N/PERTR/DQ1 ( 161 , 41 ) , DQ2(  1 61 . 41 ) , DQ3( 1 61 , 41 ) ,DQ*(161 ,41 ) 
CQNMON/GR I D 1 /X ( 1 6 1 . 41 ) , Z( 1 61 , 41 ) 

CDM40N/DGR I D/DT . IMAX . KMAX . I  TEL , I TEU 

C0**40N/FL0W/Q1(161 .41  )  .02(161 ,41  )  .03(161  ,41)  ,04(161 ,41  ) 
COLMON/PAR/GAAMA . REYREF , ALFA . ALFA  1 . REDFRE . AM I NF , ALFA I 
DIMENSION  RHS( 161 .4) 

X T AU ( I ,K)  -  OMEGA  .  Z( I ,K) 

YTAU(I.K)  -  -  OMEGA  .  X(I,K) 

THIS  SUBROUTINE  COMPUTES  THE  RESIDUAL  ON  THE  RIGHT  HAND 
SIDE  ARISING  FROM  THE  EULER-  PART  OF  THE  GOVERNING  EQUATIONS 


FLUX  TERMS  WITHIN  THE  XI-  DERIVATIVE 
DO  100  K  -  2  .  KMAX  -  1 
DO  10  I  -  1  ,  IMAX 

UCON  -  (02(1 ,K)/Q1(I.K)  )  •  ( Z ( I , K+1 )— Z( I , K— 1 ) ) 
1  -  (03(1 ,K)/01(1,K))  •  ( X( I ,K+ 1 )— X ( I , K— 1 ) ) 

UCON  -  0.25  •  DT  •  UCON 


I 


* 


v  .v  -  v 


Si 


3 


X I T  -  -  XTAU(I.K)  •  (  Z  (  I , K+ 1 ) — Z ( I . K— 1 ) ) 

1  +  YTAU(I.K)  •  ( X ( I , K+1 )  -  X( I ,K-1 ) ) 

XIT  -  X I T  .  DT  .  0.25 
UCON  -  UCON  +  XIT 
RHS(I.I)  -  UCON  •  01(1 ,K) 

R  -  1 .  /  01 (I ,K) 

P  -  ( GAMMA- 1 . )  .  (04(1, K)  -  .5  •  R  • (02 ( I ,K ) • »2+ 

1  03(1. K).. 2)) 

RHS( 1.2)  -  02(1. K)  .  UCON  +  P  .  OT  .  0.25  •  (Z(I.K+1)  -  Z(I.K-I)) 

RHS( 1.3)  -  03(1. K)  •  UCON  -  P  .  DT  •  0.25  •  ( X ( l , K+1 )-X ( I , K-1 ) ) 

RHS(  1.4)  «  UCON  .  (04(1. K)+P)  -  XIT  •  P 

10  CONTINUE 

DO  1 1  I  -  2  .  I  MAX  -  1 

DOI(I.K)  =  D0 1(1  ,  K )  -  RHS( 1  +  1.1)  +  RHS(I-I.I) 

DQ2 ( I . K )  =  DQ2 ( I . K )  -  RHS (1  +  1 .2)  +  RHS(I-1.2) 

003(1. K)  *  DQ3 ( I . K )  -  RhS( 1+1 . 3)  +  RHS(I-1.3) 

11  004(1. K)  -  004(1, K)  -  RHS( 1+1 .4)  +  RHS ( 1—1 ,4) 

100  CONTINUE 


FLUX  TERMS  WITHIN  THE  ETA-  DERIVATIVE 


DO  200  I  «  2  .  IMAX  -  1 
DO  20  K  -  1  ,  KMAX 

VCON  -  (02(1 ,K)/01 (I .K)  )  »  (2(1-1 . K )— Z ( 1+1 .K)) 

1  +(03(I.K)/Q1(I.K)  )  •  ( X ( I +1 , K  )— X( 1—1 ,  K )  ) 

VCON  =»  VCON  .  0.25  •  DT 

ETAT  -  -XTAU(I.K)  •  (Z( 1—1 ,K)-Z(I+1.K))  -  YTAU'I.K). 

1  (X( 1+1 ,K)-X(I-1 ,K)) 

ETAT  -  ETAT  .  0.25  .  DT 
VCON  -  VCON  +  ETAT 
RHS(K. 1 )  -  VCON  .  01 ( 1 ,K) 

P  -  (GAM4A-1 . )  •  (04(1. K)  -  0.5  • (02( I . K ) ••2+Q3( ].K)»»2)/01(I,K)) 
RHS(K . 2 )  -  VCON  .  02(1 . K )  +P  .  DT  .  .25  •  (Z( 1-1 , K )-Z( 1+1 . K ) ) 
RHS(K , 3)  -  VCON  «  Q3 ( I , K )  +  P  .  DT  •  . 25  •  (X(]+1.K)  -  X(I-I.K)) 
RHS (K . 4)  -  VCON  .  (Q4(I,K)+P)  -  ETAT  •  P 


RHS (K . 4) 

20  CONTINUE 
DO  21  K 
DO 1 ( I . K ) 
002(1 .K) 
DQ3( I ,K) 

21  D04(I.K) 
200  CONTINUE 

RETURN 

END 


2  KMAX  —  1 

*  DQI(I.K)  -  RHS ( K+1 . 1 )  +  RHS (K-1 . 1 ) 
»  D02 ( 1 . K )  -  RHS(K+1.2)  +  RHS(K-1 .2) 

-  DQ3( l  ,K)  -  RHS(K+1 .3)  +  RHS(K-1 ,3) 

-  D04 (  I  .  K )  -  Rl  ;S(K+1 . 4)  +  RHS(K-1.4) 


SUSROUT INE  ROTGRID(X,Z. IMAX . KMAX .DALFA) 

ROTATE  GRID  IN  THE  CLOCKWISE  DIRECTION  BY  AN  AMOUNT  DALFA 
DIMENSION  X( 1 61  .41 ) .Z( 161  .41 ) 

CA  =  COS (DALFA) 

SA  «=  -  S I N (  DALFA) 

DO  20  K  -  1  ,  KMAX 
DO  20  I  -  1  ,  IMAX 
XI  «  X( I .K) 

Z1  =  Zl I . K  ) 

X( I  ,K)  «  XI  •  CA  -  Zl  •  SA 


X( I  ,K)  «  XI  •  CA  -  Zl  •  SA 
20  Z( I . K)  -  Zl  •  CA  +  XI  •  SA 
RETURN 


SUBROUTINE  CPPL0T(I1  . I  2 . FMACH , X . Y .CP) 

THIS  SUBROUTINE  PLOTS  CP  AT  EQUAL  INTERVALS  IN  THE  MAPPED  PLANE 
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C0M40N/SK!NCF/CF(161 ) 

DIMENSION  KODE(4).LINE(90).X(161),Y(161),CP(161) 

DIMENSION  CFX(3.49).CFY(3,49),CPX(3.49),CPY(3.49) 

DATA  kOOE/IH  , 1H+, 1HI , 1H»/ 

WRITE  (6.2) 

2  FORUAT ( 50H0PLOT  OF  CP  AT  EQUAL  INTERVALS  IN  THE  MAPPED  PLANE/ 
1  10H0  X/C  , 1 0H  CPU  ,10H  CFU  ,10H  CFL 

CP0  -  (1.  +  .2  •  FMACH  »«2)  ••  3.5  -  1. 

CP0  -  CP0  /  (  .7  •  FMACH  »»2) 

K0  -  30.  •  CP0  +4.5 
IMIN  -  (I2-I1)/?  +  II 
I  LOW  -  2  •  IMIf. 

I  COUNT  -  0 
CHD-X(II)  -  X(IMIN) 

DO  12  I  -  1  .90 
12  LINE(I)  -  KODE(I) 

DO  34  1  -  IMIN  .  12 

K  -  30.  •  (CP0  -  CP(I))  +  4.5 

K 1  -  30.  •  (CP0  —CP ( I  LOW- 1 ) )  +4.5 

IF(K.LT.I)  K  -  1 

IF(KI.LT.I)  K 1  -  1 

I  F( K . GT . 90)  K  -  90 

I  F(  K 1  .GT.  90)  K 1  -  90 

LINE(K0)  -  KODE(3) 

LINE(K)  -  KODE(2) 

LINE(KI)  -  KOOE(4) 

XOC  -  (X(I)  -  X( IMIN) )  /  CHD 
WRITE  (6.610)  XOC. CP(I),CF(I),CF(1L0W-I). LINE 
L1NE(K1)  -  KODE( 1 ) 

34  LINE(K)  -  KODE( 1 ) 

..  GENERATE  PLOT3D  CP  AND  CF  PLOTTING  FILES 
DO  500  I-IMIN. 12 

XOC  -  (X(I)  -  X( IMIN) )/CHD 
I COUNT  -  I COUNT  +  1 
CPY( 1 . I COUNT)  «  0.000000 
CFY( 1 . ICOUNT)  -  0.000000 
CPY(2. ICOUNT)  -  CP(ILOW  -  I) 

CFY(2. ICOUNT)  -  CF(ILOW-I) 

CPY(3 , ICOUNT)  -  CP( I ) 

CFY(3. ICOUNT)  -  CF ( I ) 

CPX(1 . ICOUNT)  -  XOC 
CFX( 1 , ICOUNT)  -  XOC 
CPX(2. ICOUNT)  -  XOC 
CFX(2, ICOUNT)  -  XOC 
CPX(3. ICOUNT)  -  XOC 
CFX( 3. ICOUNT)  -  XOC 

0  CONTINUE 
I  DM  -  3 

JDM  «  12  -  IMIN  +  1 
WR I TE( 50)  IOM.JDM 

WR I TE( 50) ( (CPX( I , J ) .  1-1 ,IDM) .J-1 ,JDM). 

+  ((CPY(I.j).  1-1 . IDM) . J-1 . JDM) 

WR I T£ ( 60)  IDM. JDM 

WR I TE ( 60 ) ( (CFX ( I . J) .  1-1 . IDM) .J-1 .JDM)  , 

+  ((CFY(I.J),  1-1 . IDM). J-1 , JDM) 

RETURN 

610  FORMAT ( 4F 1 0 . 4 , 90 A 1 ) 


®  (9 


TITLE: 

NACA  0012  AIRFOIL 


IMAX: 

KMAX  : 

DT: 

WW: 

ALFA: 

ALFA1 : 

ALFA! 

161 

41 

.005 

5. 

15.00 

10.00 

19.00 

ITEL: 

I  TEU : 

REYREF 

:  DNMIN: 

TSTART  : 

FORMAT : 

RSTRT 

31 

127 

3  45 

00005 

-1  .0 

3.0 

TRUE 

CSTP : 

CPLT : 

NSTP  : 

PSTP : 

0 

0 

1000 

50 

FNU: 

FNL 

FSYVt: 

33.  33.  1. 

X :  Y: 

0. 

.0050  .01153 

.0125  .01894 

.0250  02615 

0500  .03555 

.0750  .04200 

.1000  .04683 

.1500  .05345 

.2000  .05737 

.2500  .05941 

.2600  .05962 

.2700  .05978 

.2800  .05992 

.2900  .05999 

.3000  .06000 

.3100  .05999 

.3200  .05992 

.3300  .05980 

.3400  05965 

.3500  .05947 

4000  .05803 

.4500  .05581 

5000  .05294 

.5500  04952 

.6000  .04563 

.6500  .04137 

.7000  .03664 

.7500  .03161 

.8000  .02623 

.8500  .02055 

.9000  01448 

.9500  00B07 

1.0000  .00126 


REDFRE:  AMINF 
0.151  .283 

PITCH:  PLUNGE: 
TRUE  FALSE 
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